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1  INTRODUCTION 

Interplanetary  missions  will  need  rockets  with  high  specific  impulse  and  high  specific  power. 
Nucleeir  fusion,  the  power  source  of  the  Sun,  could  power  a  rocket  for  these  interplanetary 
missions. [1]- [8]  While  a  voyage  to  Mars,  for  example,  would  take  almost  a  year  with  chemical 
propulsion,  a  fusion-powered  rocket  could  get  there  in  months.  The  high  temperature  of 
fusion  plasmas  (already  found  in  experimental  fusion  reactor  concepts)  would  give  /^p’s  of 
lOOO’s  of  seconds.  The  large  energy  release  from  fusion  reactions  would  give  the  necessary 
specific  power. 

If  the  rocket  could  pick  up  more  propellant  when  it  arrived  at  Mars,  it  would  not  need 
to  carry  out  propellant  for  the  voyage  back.  This  would  reduce  the  overall  mass  of  the  in¬ 
flight  rocket,  increase  the  thrust-to-weight  ratio,  and  reduce  transit  time.  A  source  of  such 
propellant  is  the  Martian  atmosphere,  which  is  95%  CO2. 

The  JANAF  Thermodynamic  Tables  [9]  give,  up  to  6000  K,  the  enthalpy,  entropy,  and 
heat  capacity  of  several  elements  and  molecules.  However,  fusion  plasmas  have  temperatures 
in  the  lOO’s  of  keV  (10®  K).  By  rederiving  theoretically,  as  the  JANAF  tables  do,  the  ther¬ 
modynamic  properties  of  the  fusion  rocket  propellant,  but  using  extended  excitation  levels, 
we  can  estimate  rocket  performance  above  6000  K. 

At  high  temperatures,  as  molecules  in  the  gas  dissociate  into  their  atomic  components, 
chemistry  simplifies.  Although  probability  increases  with  temperature  that  a  molecule  will 
occupy  a  higher  excited  state,  the  number  of  molecules  decreases  with  temperature.  Once 
most  of  the  molecules  have  dissociated,  the  remaining  molecules  have  little  effect  on  the 
thermodynamic  properties  of  the  total  gas. 

At  temiicraturcs  below  6000  K,  the  number  of  chemicals  to  consider  is  overwhelming. 
Computer  codes  exist,  such  as  Selph’s  “ISP  Code”  {10],  that  first  identify  chemicals  present 
at  a  given  temperature  and  pressure,  and  then  interpolate  JANAF  table  data  to  find  their 
partial  densities  at  equilibrium.  From  this  they  calculate  the  total  gas  thermodynamic 
properties. 

The  code  described  here  considers:  mixtures  of  carbon  (C),  hydrogen  (H),  nitrogen  (N), 
and  oxygen  (0),  up  through  H//,  Cv,  N/v,  and  Ov^;  any  diatomic  combinations  of  these 
elements,  both  neutral  and  singly  ionized.  The  user  can  enter  in  spectroscopic  data  for 
his  own  elements.  Program  results  agree  with  Selph’s  code  at  6000  K,  for  dominant  species; 
output  graphs  of  nitrogen  species  densities  (3000  K  to  30,000  K)  and  air  si)e(  ies  mole  fractions 
(3000  K  to  10,000  K)  match  Cambel[ll]. 
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2  THEORY 

2.1  Saha  Equation 

Mass  Action.  Assume  diatomic  molecules,  A2,  which  dissociate  into  atoms,  .4.  At  equi¬ 
librium,  write  the  dissociation  reaction 


A2  ^  2A  —  uo,Ai  > 


(1) 


where  represents  the  dissociation  energy  of  the  reaction.  According  to  the  law  of  mass 
action  [12],  write  the  equilibrium  constant  of  the  reaction  in  terms  of  the  inventories,  N,  of 
the  individual  species. 


K  - 


Na 


i‘^) 


Partition  Function.  The  total  energy  of  a  particle  consists  of  its  translational  and 
internal  energies.  Neglect  nuclear  excitation,  and  assume  separable  rotational,  vibrational, 
and  electronic  energies,  to  write  the  particle  Hamiltonian,  H, 

"Hiot  —  "Htr  +  Tirol  +  Tivib  +  Tig/.  (3) 

For  a  particle  partition  function,  write[13] 

Ziot  —  ZtrZrotZvibZeh  (4) 

where 


Ztr 

(5) 

Ze, 

=  Es.“p(^). 

(G) 

T 

Zrot 

Rtf  - , 

©rol 

(7) 

Zvib 

T 

(8) 

and  m  represents  the  particle  mass;  h  represents  Planck’s  constant;  k  represents  Boltzmann’s 
constant;  T  represents  the  gas  temperature;  V  represents  the  gas  volume;  u,  represents 
the  energy  of  the  ith  electronic  excitation  state;  gi  =  27-1-1  represents  the  degeneracy 
(multiplicity)  of  the  ith  electronic  excitation  state;  J  represents  the  total  angular  momentum 
quantum  number;  0^1  represents  the  characteristic  temperature  for  rotation  of  the  ground 
state  (electronic)  of  the  molecule,  and 


©rol  = 


(9) 


2 


where  I  represents  the  moment  of  inertia  of  the  molecule;  ©„,-6  represents  the  characteristic 
temperature  for  vibration  of  the  ground  state  (electronic)  of  the  molecule,  and 

0V.5  =  (10) 

where  i/  represents  the  vibration  frequency  of  the  molecule. 

Electronic  excitation  energies,  u,-,  and  the  degeneracies,  (7,,  come  from  spectroscopic 
data[14].  Similarly,  find  Qrot  and  0„jj,  for  diatomic  molecules  in  Huber  &  Herzberg(15]. 
Atomic  species  have  no  vibrational  or  rotational  energy,  so  that 


Zrot  =  Zvib  =  1  (atomic  species).  (11) 

In  terms  of  the  total  partition  function,  write  the  reaction  equilibrium  constant  in  Eq.(2) 


Ideal  Gas.  Assume  an  ideal  gas  (a  better  approximation  as  more  molecules  dissociate), 

PV  =  NkT,  (13) 

and  substitute  V  into  Eq.(5); 

^  \  -  1  - 


-'Ir 


\  /i2 ; 


■N, 


(14) 


where  N  represents  the  total  number  of  particles  in  the  gas,  specifically  “heavy”  particles, 
Nh,  (molecules,  atoms,  ions)  and  electrons,  Ne- 


N  =  Nh  +  Ne. 


1,15) 


Normalize  all  species  particle  inventories  to  the  total  number  of  nuclei,  Nq,  present  in  the 
gas,  and  define  the  nuclear  fractions  (similar  to  Cambel[ll]): 


Va,  = 
Va  = 
Vh  = 
Vc  = 


No  ' 
No' 

No' 


Using  Eqs.(16)-(19)  in  Eqs.(2)  and  (12)  gives  the  first  Saha  equation: 

3/2 


Va  _  f  ^  (fcT)^/^  Zel,A 

yAiiVh  +  ye)  V  niAj  J  P  {ZciZ  rot^vib^Aj 


exp 


f-UD,AA 
\  kT  )' 


(16) 

(17) 

(18) 

(19) 

(20) 


3 


using  the  component  partition  functions  for  the  total  partition  functions,  and  simplifying 
the  translational  partition  function  for  and  A.  Note  that  by  taking 
the  second  right-hand-side  term  (in  parentheses)  becomes  the  reduced  mass  of  the  particles. 
Also  note  that,  for  a  given  temperature  and  pressure,  the  value  of  the  right-hand-side  of 
Eq.(20)  is  known  (abbreviated  R). 


Ionization.  Now  consider  ionization  as  well  as  dissociation.  At  ctiuililjiium, 

A+-'  ^  -He"  -  7i/j, 


(21) 


where  j  +  \  represents  the  degree  of  ionization  and  u/j  represents  the  energy  of  tlie  yth 
ionization.  Equate  the  two  forms  of  the  reaction  equilibrium  constant  to  write 


VA+i+^Ve  _  /-77;,A 

yA-yiiyh  +  ye)  V  j  p  kT  )■ 


(22) 


Chemical  Reactions.  Now  consider  a  chemical  reaction  between  elements  A  and  B  that 
forms  the  compound  AB.  At  equilibrium, 


AB  ^  A  +  B  —  UR^ABy 


and  so 


yAB{yk 

If  the  compound  AB  ionizes, 


and 


AB  ^  AB'^  +  e  -  u/^ab, 
VAB+y^  _  /'27rme\^^^  (kT)^^^  {ZelZrotZ^ib)AB* 


yABiVh  +  ye)  \ 


P  {ZelZrot^vib)AB 


exp 


(  ~'^kAB\ 

\  kT  ) 


Saha  Equations.  Solve  the  following  set  of  equations: 


(23) 


yAya  _  /27ry/^  /  m^mfl  \3/2  (fcy)5/2  Zel,AZel,D  ,  (-Un^AD^  ^24^ 

liVh  +  y.)  \hV  \m^  +  mo)  P  {Z„Z,A,,)Aa\  kT  )'  '  ’ 


(25) 

(26) 


yA+ye 

=  i^Aly 

(27) 

yAivh  +  ye) 

yA-^^Ve 

=  Ra2, 

(28) 

yA^yh  +  ye) 

VA+^ye 

II 

(29) 

VA^^iVh  +  ye) 

yA+*ye 

=  Ra4, 

(30) 

yA+»{yh  +  ye) 

BA 

yA,(y/.  +  y*) 

= 

^A6. 

(31) 

Rai, 

(32) 

yAiiyn  +  ye) 

ya+ye 

Rbu 

(33) 

yaiyn  +  ye) 

ya-k^ye 

Rd2, 

(34) 

yfl+(y/>  +  ye) 

ya+^ye 

RoZy 

(35) 

yfl+s(yh  +  ye) 

yB-k*ye 

P-BAy 

(36) 

yo+i(y/.  +  ye) 

yl 

Rady 

(37) 

yflj(y/i  +  ye) 

(38) 

"2 

Raiy 

yBAyh  +  ve) 

VaVb 

Rabu 

(39) 

yAB^yh  +  ye) 

yAB+ 

RAB2y 

(40) 

yABivh  +  ye) 

yA  +  yA+  +  yA+»  +  yA+»  +  yA+«  +  y^?  +  y/t+  + 
yfl  +  yB+  +  yi?+*  +  ys**  +  yB+*  +  yaj  +  ye*  + 

(41) 

yAB  +  yAi?+ 1 

yA+  +  2y>i+i  +  3y4+3  A  4y^+4  -f  y^+  + 

3 

(42) 

yB+  +  2yij+2  +  3yB+3  +  4yg+«  +  y£j+  + 

VAB+i 


yo,A  =  yA  +  yA+  +  yA+>  +  yA+’  +  va**  +  2yAj  +  2y^+  +  uab  +  yAB+ ,  (43) 

yo,B  =  yB  +  yB*  +  yo+*  +  yB+»  +  ya**  +  2y5,  +  2y^+  +  y^£?  +  Vab-^  )  (44) 

1  =  yo,A  +  yo.oj  (45) 


where  the  R  (for  “RHS” ,  right-hand-side)  represent  the  known  functions  of  temperature  and 
pressure,  and  yo,/»  and  yo,^  represent  the  stoichiometric  fractions  of  elements  A  and  B. 
With  solved  Eqs.(27)-(45),  the  mole  fraction  of  the  fcth  species  becomes 


Nk  _  Npijk  _  Vk 

N  iVo(y/,  +  ye)  ”  (y/.  -t-  ye) ' 

Similarly,  the  partial  density  of  the  A:th  species  becomes 

,  P  Vk  P 

V  "  Aro(y/.  +  y.)kT  ■  (y,  -h  ye)  kT' 


(46) 


(47) 
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2.2  Thermodynamic  Functions 

Sum  the  energies  in  translation,  excitation,  vibration,  electronic  excitation,  dissociation, 
ionization,  and  chemical  reactions  to  get  the  total  energy,  U,  of  the  gas.  Thus, 

U  =  Utr  +  Urot  +  +  Uel  +  Ud  +  Ui  +  Ur.  (48) 


Sum  the  internal  energy  of  the  gas  and  the  thermal  excitation  to  get  the  enthalpy  of  the  gas, 
H: 


H  =  U  +  PV. 


(49) 


Divide  H  hy  n  =  N/Nav,  the  number  of  moles,  where  Nav  is  Avagadro’s  niiinbcr,  to  get  the 
molar  enthalpy,  h: 


h  =  (50) 

substiting  PV  =  NkT  into  the  second  term.  Divide  the  gas  energy  by  the  number  of  particles 
to  get  the  average  particle  energy, 

■ —  '^vib  d”  "^rot  T  d”  U/j  "h  'Uj  4"  U/j. 


Examine  each  of  the  average  energies  separately. 

Each  of  the  three  degrees  of  freedom  of  translation  has  kT/2,  and  so  the  average  trans¬ 
lational  energy,  Utr,  becomes 

u,r  5=  hr.  (52) 

Assume  that  the  average  rotational  and  vibrational  energies  for  diatomic  molecules  have 
kT  each,  a  valid  approximation  at  high  temperatures. 

Urot  =  Uvib  =  kT.  (53) 


Find  the  average  electronic  excitation  energy,  Ugi,  from 


T,iUigiexp{-Ui/kT) 

E,-5.exp(-u,/fcT) 
E,-  UjQi  exp{-UilkT) 

Zel 


(54) 


Find  the  average  dissociation,  ionization,  and  reaction  energies: 


=  Up, 

=  u/,  and 

=  Ur. 


ud 

ui 

Ur 


6 


(55) 

(56) 

(57) 


Recognize  that  Nav^  =  R  =  8.314  J/mol,  the  gas  constant,  so  the  enthalpy  becomes 


h  =  RT 


kT  )  ■ 


(58) 


Atomic  species  do  not  have  the  inner  terms  in  parentheses,  since  they  have  neither  rotational 
nor  vibration  energies,  nor  have  reacted  chemically  to  form  a  compound. 

Define  the  constant  pressure  heat  capacity,  Cp, 


Cjj  — 


±h 

dT^' 


(59) 


Recall  the  definition  of  n*/  from  Eq.(55),  to  find  the  rightmost  term; 

f  E,-  UiQi expi-Ui/kT)\ 

dT  V  z^i  ;  ’ 

=  exp(-tii/fcr)  - 

Ei  U.gi  exp(-u,/A:r)  d  ^ 

Zl  dT 

1  (J^iu‘jgiexp{-Ui/kT) 
kT^  [  Ze/ 

Ei  Ujgi  exp(-Ui/kT)  E,-  UjQj  expj-Uj/kT) ) 

Zei  Zel  } 

= 

Therefore,  write  the  heat  capacity 


(60) 


Express  a  difference  in  entropy  by[l6] 

^  PI  ^2 

Operate  on  Cp,  in  Eq.(61),  to  evaluate  the  first  right-hand-side  term: 


lnp  +  /V4u 


a  (IT 


(61) 


(62) 


(63) 
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Substitute  all  the  tenns  back  into  the  original  equation  to  get 

»2-^.=H[(|  +  (2))ln|  +  ^  +  ln|!i-ln^].  (67, 

For  convenience,  take  state  1  at  standard  temperature  and  pressure  (STP:  T  =  298.15  K, 
P  =  1  atm). 


Gas  Thermodynamic  Properties.  Equations  (58), (61),  and  (67)  give  formulae  for  the 
enthalpy,  heat  capacity,  and  entropy  difference  for  individual  species  present  in  the  gas.  To 
get  the  entropy,  heat  capacity,  and  entropy  for  the  entire  gas,  sum  the  contributions  from 
each  species.  For  example,  write  the  gas  enthalpy 

where  hi  represents  the  enthalpy  of  the  tth  species,  and  yi/{yh  +  ye)  represents  the  mole 
fraction  of  the  tth  species.  Find  gas  heat  capacity  and  enthalpy  difference  in  a  similar 
manner. 
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2.3  Rocket  Applications 

To  estimate  rocket  performance,  make  several  assumptions  [18]: 

1.  Homogeneos  propellant  does  not  vary  in  composition  throughout  the  rocket  chamber 
and  nozzle. 

2.  The  propellant  obeys  the  perfect  gas  laws. 

3.  The  propellant  flows  adiabatically. 

4.  The  propellant  flows  steadily;  expansion  occurs  iseiitropically,  without  shocks. 

5.  Exhaust  gases  have  axially  directed  velocity. 

6.  Nozzle  expansion  does  not  change  the  chemical  equilibrium  established  in  the  chamber 
(“frozen  flow”). 

7.  Propellant  contains  gaseous  species  only. 

The  diatomic  species  do  not  exacly  follow  the  ideal  gas  assumptions,  since  they  have 
degrees  of  freedom  in  rotation  and  vibration.  At  high  temperatures,  as  the  diatomic  species 
dissociate,  the  ideal  gas  assumption  improves.  As  the  nozzle  expands  and  the  temperature 
and  pressure  decrease,  the  dissociation/ionization  equilibrium  changes.  However,  assume 
that  expansion  occurs  quickly  enough  that  the  equilibrium  at  the  exit  remains  unchanged 
from  the  chamber,  and  that  the  expansion  occurs  isentropically  (reversibly).  Use  the  frozen 
flow  approximation,  rather  than  calculate  the  equilibrium  conditions  at  both  the  chamber 
and  exit,  since  the  Saha  equations  do  not  give  correct  answers  at  very  low  pressures.^ 

Since  the  expansion  occurs  isentropically, 


where  7  =  5/3  =  Cp/c„  for  an  ideal  gas.  The  heat  capacity  for  an  ideal  gas  remains  constant 
at  Cp  =  Rbf  2,  so  conservation  of  energy  yields 

(70) 

=  Cp(T„-T). 


Assume  that  the  axial  velocity  of  the  propellant  in  the  chamber,  v,  is  small  compared  to 
the  exit  velocity,  t;***  to  solve  for  the  exit  velocity: 


^The  Saha  equations  describe  equilibrium  brought  about  from  collisions.  At  low  pressures,  fewer  collisions 
occur,  and  equilbrium  occurs  by  way  of  longer-range  interactions  (“coronal  equilibrium").  This  project 
considered  only  Saha  equilibrium. 
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Convert  R  from  J/K/mol  to  J/K/kg  by  dividing  R  by  the  average  molecular  weight,  MW, 
where 

_  yo,AMWA  +  Vo.bMWb 


(Vh  +  J/e) 


to  get 


\  ^  ~  (~P“ ) 


The  exit  velocity  gives  thrust,  F,  specific  impulse,  Itp,  and  jet  power,  Pjet'. 


F 

Itp 

Piet 


mVe^, 

— ,  and 
9 


(74) 

(7o) 

(7G) 
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3  USER’S  MANUAL 


This  section  tells  how  to  use  the  code,  what  the  input  parameters  mean,  what  techniques  work 
best  for  the  desired  results,  where  to  look  for  output,  and  how  to  create  new  input  elements.  If 
the  user  has  trouble,  the  author  receives  internet  email  at  the  address  nachtriebfiuiuc .  edu, 
or  US  mail  care  of  the  Phillips  Laboratory. 

3.1  Starting  Up 

By  making  a  directory  on  the  hard-drive  to  which  all  the  program  files  may  be  copied,  the 
user  can  keep  input  and  output  files  separate  from  other  programs,  and  can  keep  track  of 
the  output  files  the  program  creates.  To  do  so,  at  the  prompt  (something  like  C:\>),  type: 

C:\>  md  saha 

C:\>  copy  a:*.*  c:\saha 

C:\>  cd  \saha 

C:\SAHA> 

These  commands  make  a  directory  named  SAHA,  copy  the  program  files  from  floppy  to  the 
new  directory,  and  then  change  into  the  new  directory. 

To  run  the  code,  type  the  program  name,  ab: 

C:\SAHA>  ab 

Several  examples  might  prove  helpful  for  demonstrating  how  to  operate  the  code. 

3.2  Example:  Carbon  Dioxide,  One  Temperature 

A  few  lines  of  introduction  greet  the  user,  after  which  the  program  prompts  the  user  to 
indicate  which  two  elements  he  wishes  the  program  to  work  with. 

Program;  AB.FOR,  ver.  1.0;  Author:  R.Nachtrieb,  Aug  92; 

OLAC-PL/RKFE  Edwards  AFB  CA  93523-5000 
internet  email:  nachtriebiOuiuc.edu 


C,  H.  N,  0 

Which  two  (2)  elements  to  use? 

CC  .0  ] 

The  user  types  the  letters  corresponding  to  the  symbols  of  the  elements  he  wishes  to 
consider.  The  program  suggests  four  elements:  carbon,  hydrogen,  nitrogen,  and  oxygen, 
as  indicated  by  C,  H,  N,  0.  In  hard  brackets,  the  program  displays  the  current  (default) 
setting.  If  the  user  wants  to  use  the  default  settings,  he  simply  presses  <Enter>.  Otherwise, 
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he  can  type  in  two  element  symbols,  separated  by  either  a  space  or  a  comma.  Assume  that 
the  user  accepts  the  default  elements  and  presses  <Enter>. 

Once  the  user  has  entered  in  his  selections,  the  program  echoes  back  the  selection.  Note 
that  the  program  is  case  sensitive,  so,  for  example,  it  does  not  recognize  n  as  nitrogen;  the 
user  should  use  N  for  nitrogen. 

Using  C  ,0  ,C0 

Stoichiometric  Proportions  of  C  ,0  ? 

C  l.OOOE+00,  2.000E+00] 

The  program  then  asks  the  user  for  the  relative  proportions  of  each  element.  For  example, 
if  the  user  wished  to  analyze  carbon  dioxide,  the  user  could  indic.ate  that  there  was  oik;  part 
of  carbon  for  two  parts  of  oxygen.  Note  that  this  prograni  doi's  not.  consider  polyatomic 
molecules,  only  diatomic;  but  the  user  can  approximate  gases  like  CO2  by  giving  the  same 
stoichiometric  proportions  as  CO2. 

Proportions  of  C  ,0  l.OOOE+00  2.000E+00 

Chamber,  exit  pressures  (atm)? 

[  l.OOOE+00,  O.OOOE+00] 

The  program  prompts  the  user  for  the  chamber  pressure,  in  atmospheres,  to  use  in  the 
Saha  equations.  The  exit  pressure  is  used  to  calculate  ideal  rocket  performance.  At  low 
pressures,  the  Saha  equations  no  longer  give  the  correct  equilibrium  species  distribution, 
since  collisions  between  gas  molecules  occur  infrequently.  Thus,  while  the  j)rogram  may  give 
answers  for  chamber  pressures  below  10"^  atm,  the  user  should  not  trust  them.  The  program 
does  not  invoke  the  Saha  equations  to  calculate  equilibrium  conditions  at  the  nozzle  exit,  so 
the  user  may  enter  zero  pressure  (vacuum)  for  the  exit. 

pressl,press2  (atm)  l.OOOE+00  O.OOOE-i-00 

Chamber  temperature  range:  tempi (K) ,temp2(K) .howmany? 

[  6 . OOOE+03 .  3 . OOOE+04 ,  1] 

The  program  asks  the  user  for  the  temperature  range  to  consider.  If  the  user  wishes  to 
consider  a  single  temperature,  he  types  in  the  low  and  high  temperatures  of  the  temperature 
range,  and  then  types  a  1.  The  program  checks  to  see  how  many  temperatures  the  user 
wants  to  consider:  if  the  user  has  entered  in  a  1,  the  program  ignores  the  second  (high) 
temperature  and  considers  only  the  first.  Thus,  to  consider  a  gas  with  temperature  GOOD  K, 
the  user  types  6000,  followed  by  a  space,  another  temperature  (any  number  will  do,  since 
the  program  ignores  it),  and  then  a  1. 

The  program  can  run  for  temperatures  below  6000  K.  However,  since  it  does  not  consider 
polyatomic  molecules,  the  program  looses  accuracy  the  further  below  6000  K  the  user  specifies 
the  temperature.  Below  6000  K,  the  user  should  use  Selph’s  “Isp”  code. 


12 


tempi (K) ,temp2(K) .howmany  6.000E+03  3.000E+04  1 

View  loop  progress:  inner, middle, outer? 

[  20,  20,  20] 

The  numbers  that  the  user  enters  tell  the  program  how  many  loop  iterations  the  program 
should  count  before  printing  a  progress  message  to  the  screen.  On  personal  computers,  the 
Saha  equations  may  take  a  long  time  to  converge,  especially  if  the  computer  does  not  have 
a  numeric  co-processor.  To  reassure  the  user  that  the  program  runs  properly,  the  program 
prints  convergance  factors  to  the  screen,  with  the  number  of  times  the  progreim  loops  has 
iterated.  Using  fast  computers,  printing  out  to  the  screen  every  loop  slows  down  the  program; 
but  using  slow  computers,  the  user  can’t  tell  if  the  program  works  or  not,  since  the  program 
takes  a  long  time  to  show  signs  of  activity. 

Sugestions.  The  first  time  the  user  runs  the  program,  try  entering  20  20  20.  If  the 
program  seems  to  iterate  the  numbers  as  fast  as  it  can,  that  is,  the  program  does  not  pause 
between  progress  updates,  try  increasing  the  numbers.  The  first  two  numbers  need  not  be 
greater  than  100,  since  the  loops  converge  in  fewer  than  100  interations. 

The  program  is  “done”  when  the  final  residual  (sum  of  all  the  absolute  values  of  the 
differences  between  the  left-hand-sides  and  the  right-hand-sides  of  the  Saha  equations)  drops 
below  €,  where  c  =  10“^.  When  the  outer  loop  residual  gets  down  around  10“®,  it  starts 
bouncing  around.  If  the  residual  doesn’t  drop  below  e  within  a  specified  number  of  outer 
loop  iterations  (the  last  of  the  three  numbers  the  user  types),  the  program  stops.  Usually 
20  outer  loop  iterations  suffices. 

Note  that  only  the  final  of  the  three  numbers  affects  the  accuracy  of  the  program  results. 
In  all  cases,  the  loop  residuals  must  drop  below  e.  However,  if  the  program  can’t  quite  get 
the  outer  loop  residual  below  c,  but  comes  close,  the  user  can  stop  calculations  after,  say,  20 
iterations. 

seein,seeinid,maxout  20  20  20 

mass  flow  rate  (kg/s)? 

[  l.OOOE+00] 

Finally,  the  program  asks  the  user  for  the  mass  flow  rate.  The  program  uses  the  mass 
flow  rate  to  determine  the  thrust  an  ideal  rocket  nozzle  produces. 

After  reading  the  mass  flow  rate,  the  program  has  all  the  information  it  needs  to  run  a 
single  temperature  case.  It  echoes  back  the  information,  and  then  starts  to  work.  First  it 
reads  in  spectroscopic  data  from  disk. 

mdot :  1 

Getting  data  from  file  C.O.in 
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Getting  data  from  file  C.l.in 
Getting  data  from  file  C_2.in 
Getting  data  from  file  C_3.in 
Getting  data  from  file  C_4.in 
Getting  data  from  file  0.0. in 
Getting  data  from  file  0.1. in 
Getting  data  from  file  0.2. in 
Getting  data  from  file  0.3. in 
Getting  data  from  file  0.4. in 

Then  the  program  starts  to  work  on  solving  the  Saha  equations.  It  writes  to  the  screen 
what  temperature  (in  K)  and  pressure  (now  in  pascals,  Pa)  it  considers.  It  also  writes  six 
columns  of  numbers,  which  tell  the  number  of  loop  iterations  and  the  residuals  for  the  three 
program  loops. 

temp(K),  pressl(Pa)  6.000E+03  1.013E+05 

in. loop  wyhX-yOX  mid. loop  wyh-1  out. loop  residual 

2.000E+01  2.546E-07 

When  the  program  has  finished  with  a  temperature,  it  prints  out  the  jet  power,  thrust, 
and  specific  impulse  that  an  ideal  rocket  would  produce.  The  program  also  tells  which  file 
to  examine  to  get  the  complete  information  it  calculated. 

Pjet  (W)  5.556E+06 

Thrust (W),Isp(s)  3.333E+03  3.399E+02 

Look  for  output  in  file  (CO. out  ) 

Press  <Enter>  to  exit  program. . . 

Single- Temperature  Output  File.  We  can  print  a  hardcopy  of  the  output  file,  CO .  out, 

by  typing  at  the  MS-DOS  prompt: 

C:\SAHA>  type  co.out  >  prn 

However,  after  several  runs,  this  may  amount  to  a  large  stack  of  paper.  An  easy  way  to  view 
the  program  on  the  screen  uses  the  “Browse”  function  of  Frailey’s  DCOM  program  (or  an 
ASCII  text  editor). 

The  first  few  lines  of  the  output  file  echo  back  the  users  input. 

Using  C  ,0  ,C0 

Stoichiometric  Proportions  of  C  .0 
[  3.333E-01.  6.667E-01] 

chamber,  exit  pressures  (atm) 

[  1.013E+05,  O.OOOE+00] 
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Temperature  range:  tempi (K) ,temp2(K) .howmany 
[  6.000E+03,  3.000E+04.  1] 

mass  flow  rate  (kg/s) 

[  l.OOOE+00] 

Next,  the  file  contains  the  partition  functions  for  all  species  considered  in  the  program. 
Partition  functions  for  atomic  species  appear  in  the  first  five  columns,  in  increasing  order 
of  ionization.  The  sixth  and  seventh  columns  contain  partition  functions  for  homonuclear 
diatomic  species,  also  in  order  of  increasing  charge.  The  three  rows  of  data  contain  the 
electronic  (labeled  el  at  the  far  left),  rotational  (rot),  and  vibrational  (vib)  partition  func¬ 
tions,  respectively.  (Note  that  the  rotational  and  vibrational  partition  functions  for  the 
atomic  species  have  l.OOOE+00.) 

Partition  function  data  for  the  heteronuclear,  diatomic  molecule  follow  the  elemental 
partition  function  data.  Only  the  diatomic  columns  contain  data,  in  the  same  electronic, 
rotational,  and  vibrational  row  order. 


Partition  Function:  C 


atomic. 

atomic. 

atomic. 

atomic. 

atomic, 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el  9.377E+00 

5.939E+00 

l.OOOE+00 

2.000E+00 

l.OOOE+00 

8.794E+01 

5.570E+01 

vib  l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

2.292E+03 

2.515E+03 

rot  l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

2.249E+00 

3.090E+00 

Partition  Function:  0 

atomic. 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el  8.947E+00 

4.016E+00 

8.604E+00 

5.646E+00 

l.OOOE+00 

8.005E+01 

3.594E+01 

vib  l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

2.886E+03 

2.466E+03 

rot  l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

2.640E+00 

2.190E+00 

Partition  Function:  CO 

atomic. 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el 

8.390E+01 

4.540E+01 

vib 

3.227E+01 

1.073E+02 

rot 

2.844E+00 

3.185E+03 

The  output  file  contains  the  right-hand-sides  of  the  Saha  equations,  found  in  Eqs.(27)- 
(45). 


Right-hand-sides  of  Saha  Equation 

atomic.  atomic,  atomic,  atomic,  atomic,  diatomic,  diatomic, 

neutral  +1  +2  +3  +4  neutral  +1 
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C  1.998E-07  5.122E-19 

O.OOOE+00  O.OOOE+00 

O.OOOE+00 

1.240E+00 

6.017E-08 

0  1.502E-09  5.827E-27 

CO 

O.OOOE+00  O.OOOE+00 

O.OOOE+00 

1.083E+01 

1.782E-04 

2.434E-08 

8.278E-10 

The  solution  to  the  Saha  equations  gives  the  nuclear  fractions,  that  is,  the  species  inven¬ 
tory  divided  by  the  total  number  of  nuclei. 

Nuclear  fractions.y 

atomic.  atomic. 

atomic.  atomic. 

atomic. 

diatomic , 

diatomic , 

neutral  +1 

+2  +3 

+4 

neutral 

+1 

C  1.265E-04  6.198E-07 

7.789E-21  O.OOOE+00 

O.OOOE+00 

9.868E-09 

2.914E-11 

0  3.068E-01  1.131E-05 

CO 

1.616E-27  O.OOOE+00 

O.OOOE+00 

1.331E-02 

3.332E-01 

7.944E-06 

6.767E-06 

yh. ye. yh+ye. beta  6.535E- 

■01  2.664E-05  6.535E- 

•01  4.076E- 

■05 

Dividing  the  nuclear  fractions  by  the  sum  of  the  electron  fraction  and  the  heavy  fraction, 
yh+ye,  gives  the  mole  fractions. 


Mole  fractions.y/ (yh+ye) 


atomic. 

atomic. 

atomic . 

atomic. 

atomic. 

diatomic. 

diatomic , 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

C 

1.935E-04 

9.484E-07 

1 . 192E-20 

O.OOOE+00 

O.OOOE+00 

1.510E-08 

4.458E-11 

0 

4.695E-01 

1.730E-05 

2.473E-27 

O.OOOE+00 

O.OOOE+00 

2.036E-02 

1.216E-05 

CO 

5.099E-01 

1.036E-05 

e- 

4.076E-05 

The  output  file  contains  the  densities  of  the  species. 

per  cubic  meter. 

Species  densities  (m"'3) 

atomic . 

atomic . 

atomic . 

atomic , 

atomic. 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

C 

2.367E+20 

1.160E+18 

1.458E+04 

O.OOOE+00 

O.OOOE+00 

1.847E+16 

5.453E+13 

0 

5.743E+23 

2.117E+19 

3.025E-03 

O.OOOE+00 

O.OOOE+00 

2.490E+22 

1.487E+19 

CO  6.236E+23  1.267E+19 

e-  4.986E+19 

One  finds  the  molar  heats  of  formation  of  the  species  by  adding  half  the  dissociation 
energy  to  the  ionization  energies,  less  the  reaction  energy.  The  output  file  lists  heats  of 
reaction  for  the  species;  subsequently,  the  enthalpies  listed  include  the  heats  of  reaction  to 
give  “absolute”  enthalpy.  The  output  file  also  lists  the  total  gas  enthalpy. 


Heat  of  formation  (J/mol) 
atomic,  atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral  +1 

+2 

+3 

+4 

neutral 

+1 

16 


c 

2.996E+05 

1.299E+06 

3.462E+06 

7.710E+06 

1.343E+07 

O.OOOE+00 

1.167E+06 

0 

2.465E+05 

1.455E+06 

4.574E+06 

9.449E+06 

1.632E+07 

O.OOOE+00 

1.158E+06 

CO 

-5.152E+05 

8.356E+05 

Enthalpies  (J/mol) 

atomic, 

atomic. 

atomic. 

atomic. 

atomic, 

diatomic, 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

C 

4.305E+05 

1.512E+06 

3.864E+06 

8.483E+06 

1.471E+07 

2.369E+05 

1.399E+06 

0 

3.745E+05 

1.686E+06 

5.080E+06 

1.038E+07 

1.785E+07 

2.310E+05 

1 . 387E+06 

CO  -2.812E+05  -2.850E+05 

e-  1.247E+05 
CO  gas  3.727E+04 


The  output  file  also  contains  species  and  gas  heat  capacities  and  entropy  differences  (from 
STP). 


Heat  capacities  (C_p)  (J/K/mol) 


atomic. 

atomic. 

atomic. 

atomic. 

atomic , 

diatomic, 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

C 

2.318E+01 

2.084E+01 

2.083E+01 

2.079E+01 

2.079E+01 

4.169E+01 

3.971E+01 

0 

2.227E+01 

2.223E+01 

2.166E+01 

2.080E+01 

2.079E+01 

4.024E+01 

4.028E+01 

CO 

4.100E+01 

4.000E+01 

e- 

2.079E+01 

co 

gas  3.219E+01 

Entropy  diff 

(from  STP) 

(J/K/mol) 

atomic. 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic , 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

C 

2.318E+01 

2.084E+01 

2.083E+01 

2.079E+01 

2.079E+01 

4.169E+01 

3.971E+01 

0 

2.227E+01 

2.223E+01 

2.166E+01 

2.080E+01 

2.079E+01 

4.024E+01 

4.028E+01 

CO 

4.100E+01 

4.000E+01 

e-  6.240E+01 
CO  gas  9.313E+01 

Finally,  the  output  file  lists  the  rocket  performance. 

Thrust (N),Isp(s)  3.333E+03  3.399E+02 
Pjet  (W)  5.556E+06 

3.3  Example:  C2irbon  Dioxide,  Multiple  Temperatures 

To  examine  multiple  temperatures,  the  user  enters  in  data  in  much  the  same  way  as  for  a 
single  temperature.  Once  again,  run  the  program. 
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C:\SAHA>  ab 

and  provide  the  program  with  the  information  it  requests: 

Program:  AB.FOR,  ver.  1.0;  Author:  R.Nachtrieb,  Aug  92; 
OLAC-PL/RKFE  Edwards  AFB  CA  93523-5000 
internet  email:  nachtriebQuiuc.edu 


C,  H,  N.  0 

Which  two  (2)  elements  to  use? 

CC  .0  ] 

Using  C  .0  ,CQ 

Stoichiometric  Proportions  of  C  ,0  ? 

[  l.OOOE+OO,  2.000E+00] 

Proportions  of  C  ,0  l.OOOE+OO  2.000E+00 

Chamber,  exit  pressures  (atm)? 

C  l.OOOE+OO.  O.OOOE+00] 

pressl,press2  (atm)  l.OOOE+OO  O.OOOE+00 

Chamber  temperature  range:  tempi (K) ,temp2(K) .howmany? 

[  6.000E+03,  3.000E+04,  1] 

Now  the  user,  instead  of  selecting  only  one  temperature,  inputs  the  temperature  range 
he  wishes  to  consider,  and  the  number  of  divisions  within  that  temperature  range.  If  the 
user  enters  in  2  for  howmany  divisions,  the  program  considers  the  lower  and  upper  bounds; 
if  the  user  enters  3,  the  program  considers  the  lower  and  upper  bounds,  and  a  tcnipcratmo 
midway  between  the  two,  and  so  on. 

Try  entering  in 

6e3  30e3  3 

The  program  echoes  back  the  user’s  selection,  and  then  prompts  for  a  new  input: 
tempi (K) ,temp2(K) .howmany  6.000E+03  3.000E+04  3 

Print  density  (1)  or  mole  fraction  (2)  ? 

[1] 

When  the  user  selects  more  than  one  temperature  for  the  program  to  consider,  the  program 
outputs  information  in  a  different  way:  instead  of  writing  all  information  to  one  file,  the 
program  write  specific  information  to  many  files.  This  makes  it  easy  for  the  user  to  import 
a  file  of,  say,  species  densities  vs.  temperature  into  a  graphing  program.  In  this  way,  the 
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user  simply  looks  at  the  particular  file  to  see  enthalpy,  heat  capacity,  entropy,  mole  fraction, 
density,  or  rocket  performance.  Unfortunately,  MS-DOS  limits  the  number  of  files  a  user 
(actually,  the  program)  can  have  open  at  one  time.  Thus  the  user  must  choose  whether  to 
make  files  that  contain  species  densities,  or  to  make  files  that  contain  species  mole  fractions. 
Note  that  if  the  user  wants  both  density  and  mole  fraction,  he  can  run  the  program  twice, 
and  specify  density  the  first  time  and  mole  fraction  the  second  time. 

The  rest  of  the  inputs  follow  the  single  temperature  example. 

densmol  [1] 


View  loop  progress;  inner, middle, outer? 

[  100,  100,  20] 

seein, seemid,maxout  100  100  20 

mass  flow  rate  (kg/s)? 

[  l.OOOE+00] 
mdot ;  1 

Afier  the  user  has  input  the  information,  the  program  reads  atomic  spectra  data  from 
input  files,  and  starts  to  work. 

Getting  data  from  file  C.O.in 
Getting  data  from  file  C.l.in 
Getting  data  from  file  C_2.in 
Getting  data  from  file  C_3.in 
Getting  data  from  file  C_4.in 
Getting  data  from  file  0_0.in 
Getting  data  from  file  C.l.in 
Getting  data  from  file  0_2.in 
Getting  data  from  file  0_3.in 
Getting  data  from  file  0_4.in 

Since  the  program  now  considers  multiple  temperatures,  it  gives,  for  each  temperature, 
the  pressure,  progress  updates  for  Uie  loop  calculations,  and  rocket  performance. 

temp(K),  pressl(Pa)  6.000E+03  1.013E+05 

in. loop  wyhX-yOX  mid. loop  wyh-1  out. loop  residual 

2.000E+01  2.546E-07 

Pjet  (W)  5.556E+06 

Thrust(N),Isp(s)  3.333E+03  3.399E+02 

^temp(K),  pressl(Pa)  1.800E+04  1.013E+05 

in. loop  wyhX-yOX  mid. loop  wyh-1  out. loop  residual 
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2.000E+01  1.126E-06 

Pjet  (W)  4.561E+07 

Thrust (N),l8p(8)  9.551E+03  9.740E+02 

teinp(K),  pressl(Pa)  3.000E+04  1.013E+05 

in. loop  wyhX-yOX  mid. loop  wyh-1  out. loop  residual 

2.000E+01  1.859E-06 

Pjet  (W)  9.779E+07 

Thrust (N),Isp(s)  1.398E+04  1.426E+03 

Finally,  the  program  tells  the  user  which  files  to  examine  for  the  output  data.  Note  that 
the  filenames  follow  mnemonic  convention:  to  see  a  file  containing  rocket  thrust,  specific 
impulse,  and  jet  power  as  a  function  of  chamber  temperatme,  .see  file  C0_rkt.out;  for 
carbon  species  densities,  see  file  C_den.out,  and  so  on. 

Look  for  output  in  files: 

(C_eth.out  ) 

(O.eth.out  ) 

(C0_eth.out  ) 

(C_Cp.out  ) 

(O.Cp.out  ) 

(CO.Cp.out  ) 

(C.etr.out  ) 

(O.etr.out  ) 

(CO.etr.out  ) 

(C.den.out  ) 

(O.den.out  ) 

(CO.den.out  ) 

(CO.rkt.out  ) 

Press  <Enter>  to  exit  program. . . 

Multiple-Temperature  Output  Files.  When  the  users  selects  multiple  temperatures, 
the  output  files  created  contain  data  in  columns.  For  example,  the  file  C0_rkt .  out  contains 
rocket  performance  data,  as  a  function  of  temperature: 

CO  rocket  performance  parameters 

temp(K)  Pjet(W)  thrustCN)  Isp(s) 

6.000E+03  5.556E+06  3.333E+03  3.399E+02 

1.800E+04  4.561E+07  9.551E+03  9.740E+02 

3.000E+04  9.779E+07  1.398E+04  1.426E+03 

The  program  creates  three  output  files  for  density  (or  mole  fraction),  enthalpy,  heat 
capacity,  and  entropy.  Thus,  file  C_den .  out  contains  the  desities  of  all  the  carbon  atomic 
species,  the  homonuclear  diatomic  carbon  species,  and  the  total  electron  density.  Similarly, 


I 
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file  O.den .  out  contains  density  for  oxygen  species,  and  the  total  electron  density.  The  file 
C0_den .  out  contains  density  of  carbon  monoxide  species,  and  the  total  electron  density.  The 
program  writes  the  same  total  electron  densities  to  each  file  as  a  safety  check:  if  the  user 
plots  all  the  densities  on  the  same  graph,  discrepancies  in  electron  density  alert  the  user  to 
an  error  (for  example,  plotting  a  file  from  a  different  program  run). 

3.4  Example:  Nitrogen,  Multiple  Temperatures 

Suppose  the  user  wants  to  consider  a  single  element,  say  nitrogen.  Run  the  program,  and 
after  it  asks  for  the  elements  to  consider, 

C.  H.  N.  0 

Which  two  (2)  elements  to  use? 

[C  .0  ] 

enter  in 

N  N 

The  program  echoes  back  the  users  selection,  and  asks  for  the  stoichiometric  proportions. 
Using  N  ,N  ,NN 

Stoichiometric  Proportions  of  N  ,N  ? 

[  l.OOOE+00,  2.000E+00] 

At  this  point,  the  user  can  enter  any  proportions  he  desires,  since  both  species  arc  the  same. 
However,  a  much  better  way  exists: 

Hint.  When  considering  only  one  element,  type  0  for  one  of  the  proportions.  This  speeds 
up  the  code  considerably,  since  it  no  longer  has  to  consider  “heteronuclear”  diatomic  species. 
Try 

1  0 

after  which  the  program  responds, 

Proportions  of  N  ,N  O.OOOE+00  l.OOOE+00 
After  this,  the  program  functions  as  before. 


21 


3.5  Example:  New  Element 

Suppose  the  user  wants  to  run  the  program  with  an  element  the  program  doesn’t  yet  recog¬ 
nize.  The  program  provides  a  way  for  the  user  to  enter  the  information  the  program  needs 
to  run,  and  save  that  information  to  the  input  files. 

The  user  runs  the  program,  and  enters  in  the  elements  he  wishes  to  consider,  C  and  Pn.^ 

C.  H.  N.  0 

Which  two  (2)  elements  to  use? 

[N  ,N  ]  C  Pn 

Oops!  Couldn't  find  Pn  in  file  (ABX.in) 

Make  a  new  element?  (l=yes) 

[1] 

The  file  ABX.in  contains  input  information  about  the  elements  the  program  recognizes,  and 
the  homonucleax  diatomic  molecules  formed  by  those  elements;  specifically,  the  program 
looks  in  ABX.in  for  the  highest  ionization  state  and  the  mass  of  the  atomic  species,  for  dis¬ 
sociation  and  ionization  energy  of  the  homonucleax  diatomic  molecule,  and  for  characteristic 
temperatures  of  rotation  and  vibration.  The  program  looks  through  the  input  file  ABX .  in 
for  the  element  Pn,  but  doesn’t  find  it.  It  then  asks  the  user  if  he  wants  to  define  a  new 
element,  and  provide  its  characteristic  information.  If  the  user  says  no  (0),  the  program  asks 
the  user  for  two  new  elements.  If  the  user  says  yes  (1),  the  program  continues. 

Let ’ s  make  a  new  element . . . 

Highest  atomic  state  (1+highest  ionization  state)? 

CO] 

The  highest  atomic  state  refers  to  the  highest  ionization  level  to  be  considered  by  the  pro¬ 
gram.  State  1  means  the  program  considers  only  the  neutral  atom;  state  2  means  the  program 
considers  the  neutral  atom  and  the  first  ionized  atom,  and  so  on.  Without  modification,  the 
program  considers  up  to  state  5,  four  ionizations. 

The  program  then  asks  for  the  element  molecular  weight  (in  amu),  the  dissociation  and 
ionization  energies  of  the  homonuclear  diatomic  molecule  (in  eV).  If  the  diatomic  molecule 
does  not  form,  simply  enter  a  zero  dissociation  energy. 

ist;  [0] 

Atomic  mass  (amu)? 

[  O.OOOE+00] 

MW:  [  O.OOOE+00] 


^Pandamonium 
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Dissociation  energy  of  diatomic,  homonuclear  molecule  (eV)  ? 

[  O.OOOE+00] 
uD:  [  O.OOOE+00] 

Ionization  energy  of  diatomic,  homonuclecir  molecule  (eV)  ? 

[  O.OOOE+00] 
ul:  [  O.OOOE+00] 

Next  the  program  asks  for  the  characteristic  temperatures  of  rotation  and  vibration  for 
the  neutral  and  singly  ionized  homonuclear  diatomic  molecule,  at  the  ground  electronic  state. 
These  can  be  found  in  Huber  &  Herzberg  [15]  for  many  diatomic  molecules.  If  no  diatomic 
molecule  exists,  the  characteristic  temperatures  must  still  be  non-zero,  since  the  program 
divides  by  them  to  get  the  rotational  and  vibrational  partition  functions. 

Characteristic  temperatures,  at  ground  electronicstate  (K)  ? 

Neutral:  Rotation,  Vibration;  Singly  Ionized:  Rotation,  Vibration 
[  l.OOOE+00  l.OOOE+00  l.OOOE+00  l.OOOE+00] 

The  program  then  summarizes  the  information  just  typed  in  about  the  element  and  asks 
if  it  meets  the  users  satisfaction.  If  not,  the  user  has  the  opportunity  to  re-enter  it;  if  so,  the 
user  can  append  the  new  element  to  the  file  ABX .  in,  so  the  element  will  be  “known”  by  the 
program  the  next  time  the  element  is  used. 

charT:  [  l.OOOE+00  l.OOOE+00  l.OOOE+00  l.OOOE+00] 

Summary : 

Highest  atomic  state  (1+highest  ionization  state) 

[0] 

Atomic  mass  (eimu) 

C  O.OOOE+00] 

Dissociation  energy  of  diatomic,  homonuclear  molecule  (eV) 

[  O.OOOE+00] 

Ionization  energy  of  diatomic,  homonuclear  molecule  (eV) 

C  O.OOOE+00] 

Characteristic  temperatures,  at  ground  electronicstate  (K) 

Neutral:  Rotation,  Vibration;  Singly  Ionized:  Rotation,  Vibration 
[  l.OOOE+00  l.OOOE+00  l.OOOE+00  l.OOOE+00] 

Is  everything  ok?  (l=ok,other=redo) 

[1] 

Good.  Append  data  to  file  (ABX. in)?  (0=no,l=yes) 

[03 
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Next,  the  program  searches  through  the  file  ABXY.in  for  the  compound  foniu^d  by  tin; 
two  chosen  elements.  File  ABXY.in  contains  energies  of  dissociation  and  ionization,  and 
characteristic  temperatures  for  rotation  and  vibration,  for  heteronuclear  diatomic  molecules. 
If  the  program  cannot  find  that  compound,  it  prompts  the  user  for  the  requisite  information. 

Will  not  append  data  to  file... 

Oops!  Couldn't  find  CPn  in  file  (ABXY.in) 

Make  a  new  compound?  (l=yes) 

[1] 

Let ' s  make  a  new  compound . . . 

Reaction  energy  of  diatomic,  heteronuclear  molecule  (eV)  ? 

C  O.OOOE+00] 

RAB:  [  O.OOOE+00] 

Ionization  energy  of  diatomic,  heteronuclear  molecule  (eV)  ? 

C  O.OOOE+00] 

lAB:  C  O.OOOE+00] 

Characteristic  temperatures,  at  ground  electronicstate  (K)  ? 

Neutral:  Rotation,  Vibration;  Singly  Ionized:  Rotation,  Vibration 
[  l.OOOE+00  l.OOOE+00  l.OOOE+00  l.OOOE+00] 

As  for  the  homonuclear  diatomic  molecule,  if  the  two  elements  form  no  compound,  the 
user  enters  a  zero  reaction  energy.  The  characteristic  temperatures  for  rotation  and  vibration 
must  be  non-zero,  though. 

The  program  then  summarizes  the  information,  asks  if  it  is  satisfactory,  and  offers  to 
append  the  new  information  to  the  file  ABXY.in. 

TAB:  [  l.OOOE+00  l.OOOE+00  l.OOOE+00  l.OOOE+00] 

Reaction  energy  of  diatomic,  heteronuclear  molecule  (eV) 

[  O.OOOE+00] 

Ionization  energy  of  diatomic,  heteronuclear  molecule  (eV) 

[  O.OOOE+00] 

Characteristic  temperatures,  at  ground  electronicstate  (K) 

Neutral:  Rotation,  Vibration;  Singly  Ionized:  Rotation,  Vibration 
[  l.OOOE+00  l.OOOE+00  l.OOOE+00  l.OOOE+00] 

Is  everything  ok?  (l=ok,other=redo) 
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[1] 

Good.  Append  data  to  file  (ABXY.in)?  (0=no,l=yes) 

[0] 

Finally,  the  program  has  information  it  needs  about  the  atomie  and  molecidar  species, 
and  so  it  confines  with  the  rest  of  the  user  input. 

Will  not  append  data  to  file... 

Using  C  ,Pn,CPn 

Stoichiometric  Proportions  of  C  ,Pn? 

C  O.OOOE+00,  l.OOOE+00] 

Atomic  Spectra  Input.  The  program  reads  the  atomic  spectra  data  in  from  files  on 
disk,  which  must  adhere  to  this  naming  convention: 

filename  =  (element) _ (charge)  .  in 

where  (element)  represents  for  the  element  symbol,  and  (charge)  represents  the  charge  of 
the  atom.  For  example,  if  the  user  wanted  to  examine  Pn  up  to  the  fifth  ionization  state,  he 
would  need  to  enter  data  into  four  files,  named  Pn.O.in  through  Pn_4.  in. 

Each  file  contains  data  in  the  following  fashion:  the  first  line  contains  the  ionization 
energy  of  the  atom,  in  inverse  centimeters.  Subsequent  lines  contain  pairs  of  numbers:  the 
degeneracy  (dimensionless)  and  the  electronic  excitation  energy,  also  in  inverse  centimeters. 
The  electronic  excitation  energies  increase  up  to  the  ionization  energy  level.  The  last  line 
contains  two  negative  numbers,  to  flag  the  computer  to  stop  reading  degeneracy,  energy  level 
I^airs.  The  user  can  find  degeneracies  and  electronic  cxitation  levels  for  atomic  species  in 
Moore[14]. 

3.6  BenchmEirking  the  Code 

Internal  Checks.  The  user  can  try  to  fool  the  code  by  entering  two  identical  elements 
with  non-zero  total  nuclear  fractions.  For  example,  say  the  user  enters  nitrogen  and  nitrogen. 

Using  N  ,N  ,NN 

Stoichiometric  Proportions  of  N  ,N 
t  5.000E-01,  5.000E-01] 

The  program  does  not  recognize  that  these  two  elements  arc  different.  Call  the  two 
elements  N  and  N’.  If  the  program  works  properly,  the  program  should  predict  the  same 
thermodynamic  properties  for  elements  N  and  N’,  since  they  really  are  the  same.  .\lso,  the 
program  should  predict  equal  quantities  of  the  diatomic  molecules  N2,  N’2,  and  NN’.  The 
output  file  NN.out,  partially  listed  here,  shows  that  the  program  predicts  correctly. 
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Partition  Function:  N 


atomic. 

atomic. 

atomic , 

atomic. 

atomic, 

diatomic , 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

el  4.105E+00 

8.937E+00 

5.836E+00 

l.OOOE+00 

O.OOOE+00 

1.685E+01 

3.669E+01W 

vib  l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

O.OOOE+00 

2 . 087E+03 

2.159E+03 

rot  l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

O.OOOE+00 

1 . 768E+00 

1.890E+00 

Partition  Function:  N 

atomic, 

atomic. 

atomic. 

atomic. 

atomic , 

diatomic , 

diatomic , 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el  4.105E+00 

8.937E+00 

5.836E+00 

l.OOOE+00 

O.OOOE+00 

1.685E+01 

3.669E+01 

vib  l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

O.OOOE+00 

2.087E+03 

2.159E+03 

rot  l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

O.OOOE+00 

1.768E+00 

1.890E+00 

Partition  Function:  NN 

atomic. 

atomic. 

atomic , 

atomic , 

atomic , 

diatomic , 

diatomic , 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el 

1.685E+01 

3.669E+01 

vib 

2.408E+01 

1.037E+02 

rot 

2.779E+00 

3.175E+03 

Right'hand-sides  of  Saha  Equation 

atomic. 

atomic. 

atomic. 

atomic. 

atomic, 

diatomic. 

diatomic , 

neutral 

+1 

+2 

+3 

+4 

neutral 

A 

K  1.204E-09 

8.032E-23 

O.OOOE+00 

O.OOOE+00 

O.OOOE+00 

2.319E-03 

2.114E>10W 

N  1.204E-09 

8.032E-23 

O.OOOE+00 

O.OOOE+00 

O.OOOE+00 

2.319E-03 

2.114E-10 

NN 

2.319E-03 

2.114E-10 

Nuclear  fractions, y 

atomic. 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

N  1.390E-02 

1.027E-06 

5.069E-24 

O.OOOE+00 

O.OOOE+00 

1.620E-01 

2.104E-06 

N  1.390E-02 

1.027E-06 

5.069E-24 

O.OOOE+00 

O.OOOE+00 

1.620E-01 

2.104E-06 

NN 

1.620E-01 

2.104E-06 

yh,ye,yh+ye,beta  5.139E-01  8.367E-06  5.139E-01  1.628E-05 

Nitrogen  Density  Compare  the  density  vs.  temperature  profiles  of  nitrogen,  at  one 
atmosphere  pressure,  in  Figures  (1)  and  (2)[ll,  p.l38]  They  follow  each  other  closely,  except 
Cambel  does  not  consider  the  Nj  species. 

Air  Mole  Fraction  Compare  the  mole  fraction  vs.  temperature  profiles  of  air  (80% 
nitrogen,  20%  oxygen),  at  one  atmosphere  pressure,  in  Figures  (3)  and  (4)[11,  p.86]  Cambel 
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THn|«ratM«  r.*K  X  W 


Figure  2:  Equilibrium  chemical  compositiou  of  nitrogen  plasma  at  atmospheric  pressure. 
(Cambel,  p.l38,  used  without  permission) 
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considers  species  that  this  program  does  not:  N2O,  N2O,  0".  These  species  effect  the 
electron  density,  so  the  program  results  for  electron  density  do  not  match  well  with  Cambel, 
even  at  10,000  K.  The  dominant  species  at  10,000  K,  N,  O,  and  N2,  match  within  10%  or 
so.  However,  Camabel  keeps  the  density  constant  at  the  sea-level  value,  while  the  program 
keeps  the  pressure  constant  at  one  atmosphere. 


4000  6000  8000  10000 

lemperalure  (K) 

Figure  3:  Mole  fraction  equilibrium  chemical  composition  of  air  at  one  atmosphere  pressure. 


Temperature,  ’K 

Figure  4:  Mole  fraction  equilibrium  chemical  composition  of  air  at  sea-level  density.(Cambel, 
p.86,  used  without  permission) 


Carbon  Dioxide  The  user  should  use  Selph’s  “Isp”  Code  for  chamber  temperatures  below 
6000  K.  Selph’s  code  considers  all  chemical  reactions  and  all  relavant  species.  However,  the 
JANAF  tables  stop  at  6000  K,  so  at  temperatures  above  6000  K,  Sclph  holds  the  heat 
capacity  of  the  gas  constant  and  extrapolates.  The  author’s  program  calculates  theoretically 


the  equilibrium  state  function  using  the  Saha  Equation,  but  neglects  polyatomic  moIccul(;s. 
At  temperatures  much  above  6000  K,  this  assumption  becomes  minor,  as  almost  all  of  the 
molecules  have  dissociated. 

Table  (1)  compares  the  dominant  species  densities  predicted  by  Sel])h’s  and  the  author’s 
codes,  for  one  part  carbon  and  two  parts  oxygen,  at  6000  K  chamber  temperature  and 
one  atmosphere  chamber  pressure  pressure.  Note  that  the  dominant  species  are  close,  but 
discrepancies  arise  to  give  around  10%  difference  (for  example,  atomic  carbon,  diatomic 
oxygen,  and  carbon  dioxide).  Also  note  that  at  6000  K,  there  Selph  predicts  less  than  1% 
CO2. 


Table  1:  Particle  density  comparision  between  Selph’s  and  author’s  codes  for  CO2,  at  GOOD 
K  and  one  atmosphere.  _ 


Selph,  “Isp” 
(m-^) 

Author,  “AB” 
(m-“) 

0 

6.1E-I-23 

5.7E-)-23 

CO 

5.9E-I-23 

6.2E+23 

c 

l.OE-1-22 

2.3E+20 

^2 

8.0E-I-20 

2.4E-f22 

ICO, 

2.2E-1-20 

O.OE+00 

Thermodynamic  Functions  Table  2  compares  thermodynamic  properties  for  nitrogen 
species  N,  N'*',  and  N2,  using  the  JANAF  tables  and  the  author’s  code.  Results  match  closely 
for  atomic  species,  but  differ  slightly  for  diatomic  species. 


Table  2:  Thermodynamic  properties  comparison  between  JANAF  tables  and  author’s  pro¬ 
gram,  for  some  nitrogen  species. 


JANAF 

Author 

N 

Cp,T  =  3000  K  (J/mol/K) 

20.98 

20.96 

Cp,r  =  6000  K  (J/mol/K) 

25.54 

25.52 

As,  6000  K,  3000  K  (J/mol) 

15.63 

15.62 

Ah,  6000  K,  3000  K  (kJ/mol) 

68.42 

68.40 

N+ 

Cp,r  =  3000  K  (J/mol/K) 

20.97 

20.96 

Cp,r  =  6000  K  (J/mol/K) 

22.39 

22.37 

As,  6000  K,  3000  K  (J/mol) 

14.95 

14.93 

Ah,  6000  K,  3000  K  (kJ/mol) 

64.95 

65.00 

N2 

Cp,r  =  3000  K  (J/mol/K) 

37.05 

37.74 

Cp,r  =  6000  K  (J/mol/K) 

38.30 

43.42 

As,  6000  K,  3000  K  (J/moI) 

26.11 

28.38 

Ah,  6000  K,  3000  K  (kJ/mol) 

113.20 

124.20 
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^  A  FORTRAN  Source  Code 
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A.l  ab.for 

PROGRAM  AB 

IMPLICIT  REAL*8  (a-h.l-z) 

CHARACTER  chA*2.chB*2,name*4 

INTEGER  ist A . istB . seein . seemid , maxout , densmol , howmany 
REAL*8 

Apressl , press2 . temp . tempi , temp2 , delt . 
ftMWA.DA2.IA2.TA2(4).ionA(5) .guA(5,200,2) , 

AMWB,DB2.IB2.TB2(4) .ionB(5) .guB(5.200.2) . 
ftZ0A(7,3) .ZA(7,3) .ubarA(7) .uubA(7) .rhsA(7) .yOA,yA(7) . 
tZ0B(7,3) .ZB(7,3) .ubarB(7) ,uubB(7) .rhsB(7) ,yOB.yB(7) , 

AenthA(7) ,entrA(7) ,CpA(7) , 

AeiithB(7)  ,entrB(7)  ,CpB(7) , 
mB.IAB.TAB(4), 

tZ0AB(2.3) .ZAB(2.3) .ubarAB(2) ,uubAB(2) .rhsAB(2) .yAB(2) , 
ftenthAB(7) ,entrAB(7) ,CpAB(7) ,enthe,Cpe,entre,ye,yh,beta, 

Amdot , P j  et . thrust , Isp 
DATA  tmpO  /2.9815E2/ 

INCLUDE  format. for 

901  FORMAT  (5x, 'in. loop’ ,3x, ’uyhX-yOX' ,3x, 'mid. loop' ,6x, ’wyh-1’ , 

A3x, 'out. loop’ ,3x, ’residual’) 

WRITEC*,*) ’Program:  AB.FOR,  ver.  1.0;  Author:  ’// 
t’R.Nachtrieb,  Aug  92;’ 

WRITE(*,*)’OUC-PL/RKFE  Edwards  AFB  CA  93523-5000’ 

WRITEC*,*) ’internet  email:  nachtriebCuiuc.edu’ 

WRITEC*,*) 

WRITEC*. ♦) 

c - startup  function  gets  temperatures,  pressures,  spectroscopic 

c - data  from  user  and  files 

CALL  startup ( 

Apressl . press2 .tempi , temp2 , howmany , 

AchA , ist A , MVA , DA2 , I A2 , TA2 , guA , ionA , yOA , 

AchB , istB , MWB , 0B2 , IB2 , TB2 , guB , ionB , yOB , 

Aname , RAB , lAB , TAB , seein , seemid .maxout , mdot , densmol) 

c - get  partition  functions  at  STP,  for  later  use  in  entropy  difference 

CALL  Z1 (tmpO , ist A , guA , TA2 , ionA , ZOA . ubar A , uubA) 

CALL  Z1 (tmpO , istB , guB . TB2 , ionB , ZOB , ubarB , uubB) 

CALL  Z2 (tmpO , ZOA , ubar A , uubA , ZOB , ubarB , uubb , TAB , ZOAB , ubar AB , uubAB) 
delt=0 . 

IF  (howmany. NE. 1)  delt=(temp2-templ)/DFL0AT(howmany-l) 
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temp=templ 

DO  100  itemp=l,howmany 

WRITE(*,192) ’  temp(K),  pressl(Pa) ’ .temp.pressl 
WRITE(*,901) 

WRITEC*,*) 

c - get  equilibrium  molecule/atom/ion  fractions  using  partition  function 

CALL  saha( 

ft  pressl , temp . ist A , istB , seein , seemid . maxout , 
ft  MWA,DA2,IA2.TA2,ionA.guA, 
ft  MWB.DB2.IB2.TB2.ionB,guB. 
ft  ZOA,ZA,ubarA,uubA,rhsA,yOA,yA, 
ft  ZOB,ZB,ubarB,uubB,rhsB,yOB,yB, 
ft  enthA,entrA,CpA, 
ft  enthB.entrB.CpB, 

ft  RAB , TAB , TAB , ZOAB . ZAB . ubar AB , uubAB , rhs AB , yAB , 


ft  enthAB , entrAB , CpAB , enthe , Cpe , entre , ye . yh , beta) 

c - calculate  ideal  rocket  performance 

CALL  rocket ( 

ft  yOA, yOB,MWA,MWB, ye, yh, temp, pressl, press2,mdot .thrust , Isp.Pjet) 

c - write  all  accumulated  data  to  file  for  later  digestion 

CALL  output ( 


ft  yOA , yOB , tempi , temp2 , howmany , pressl , press2 , mdot , temp , it emp , 
ft  chA , DA2 , IA2 , ionA , ZA , rhsA , yA , enthA , CpA , entrA , 
ft  chB , DB2 , IB2 , ionB , ZB , rhsB , yB , enthB , CpB , entrB , 
ft  name , RAB , I AB , ZAB , rhsAB , yAB , enthAB , CpAB , entrAB .enthe , Cpe , entre , 
ft  P jet, thrust, Isp.densmol) 

WRITEC*,*) 
temp=temp+delt 
100  CONTINUE 

c - close  all  open  files,  tell  user  which  files  to  examine  for  output 

CALL  shutdown (howmany , chA , chB , name , densmol ) 

END  Imain  AB 

INCLUDE  STARTUP. FOR 
INCLUDE  GETX.FOR 
INCLUDE  GETXY.FOR 
INCLUDE  NEWX.FOR 
INCLUDE  NEWXY.FOR 
INCLUDE  FILENAME. FOR 
INCLUDE  OPUN.FOR 
INCLUDE  CLOZE. FOR 
INCLUDE  GETGU.FOR 
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INaUDE  Zl.FOR 
INCLUDE  Z2.FOR 
INCLUDE  SAHA. FOR 
INCLUDE  RHSl.FOR 
INCLUDE  RHS2.F0R 
INCLUDE  THERNl.FOR 
INCLUDE  THERM2.F0R 
INCLUDE  THERME.FOR 
INCLUDE  GESl.FOR 
INCLUDE  ITERY.FOR 
INCLUDE  GETYE.FOR 
INCLUDE  GETYH.FOR 
INCLUDE  BALANCE. FOR 
INCLUDE  NUY.FQR 
INCLUDE  GETWYH.FOR 
INCLUDE  RESl.FOR 
INCLUDE  RES2.FOR 
INCLUDE  ROCKET. FOR 
INCLUDE  OUTPUT. FOR 
INCLUDE  SPECSHOW.FOR 
INCLUDE  ZISHOW.FOR 
INCLUDE  Z2SH0W.F0R 
INCLUDE  OPUNHDR.FOR 
INCLUDE  SHUTDOWN. FOR 
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A. 2  startup.for 


Q^iihi^^iif:tfilHlf^:li:it:it(*^i*HHfiliit4:************************************************** 

SUBROUTINE  startup ( 
ftpr es  s l,press2,t  emp 1 . t emp2 , ho vmany . 
tchA , ist A , MWA , DA2 , IA2 , TA2 , guA , ionA ,yOA , 
ftchB . istB , MWB , DB2 , IB2 , TB2 , guB , ionB , yOB , 

Aname , RAB , I AB , TAB , see in , seemid , maxout . mdot , densmol ) 

CHARACTER  chA*2 , chB*2 , chdum*2 , name*4 
INTEGER 
AistA, 

AistB. 

Aseein , seemid . maxout , densmol . howmany 
REALMS  pressl,press2, tempi. temp2. 

AMWA. DA2,IA2,TA2(4) ,guA(5.200,2) ,ionA(5) ,yOA, 

AMWB .  DB2 , IB2 . TB2 (4) , guB (5 , 200 . 2) . ionB(5) , yOB , 

ARAB,IAB,TAB(4) 

CHARACTER  str*30 
IMPLICIT  REAL+8  (a-h,l-z) 

INCLUDE  format. for 
DATA  atm  /0.101325e6/ 

c - Read  in  default  inputs  (stored  from  last  run;  if  no  file  AB.in  exists 

c - i.e.,  first  time  the  program  is  run,  program  goes  to  label  100  and 

c - starts  with  the  “hard-wired”  defaults. 

c - User  can  select  default  by  pressing  <Enter>,  or  enter  his  own  selection. 

c - After  entering  data,  choices  eure  written  to  file  AB.in  for  next  time. 

OPEN  (UNIT=1,  FILE=’ AB.IN') 

READ ( 1 , * , ERR=100 , END=100) chA . chB 
READ(l,*,ERR=100,END=100)yOA,yOB 
READd ,  ♦ ,  ERR=100 ,  END=100)pressl  ,press2 
READCl , ♦ , ERR=100 , END=100) tempi ,temp2 , howmany 
READ ( 1 , ♦ , ERR=100 , END=100) densmol 
READ ( 1 , * , ERR=100 , END=100) seein , seemid , maxout 
READCl , ♦ , ERR=100 , END=100)mdot 
GOTO  200 
100  CONTINUE 
chA='C' 
chB='0’ 
y0A=l . 
yOB=2. 
pressl=l. 
press2=0. 
templ=6e3 
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te]ap2-30e3 

hovmanysl 

densmol«l 

seeins20 

see]&ids20 

mazout=20 

mdot=l . 

200  CONTINUE 


c - get  element  choices;  case  sensitive  (thinks  lover  case  is  a  new 

c - element,  and  prompts  user  for  element  data 


300  WRITE(*,^)'C,  H,  N,  0' 

WRITE(*,*) ’Which  two  (2)  elements  to  use?’ 
WRITE(*,’(”  ”)’)chA,chB 

READ(*,’(A)’)str 

IF  (str.NE.’  ’)  READ(str,*)chA.chB 
IF  (chA.GT.chB)  THEN 
chdum=chA 
chA=chB 
chB=chdum 
ENDIF 

OPEN  (UNIT=10,  FILE=’ABX.in’.STATUS=’old') 


c - look  in  file  ABX.in  for  information  about  generic  unknown  element 

c - “X”,  where  X  is  letter  user  typed.  If  cannot  find  “X”,  in  file 

c - ABX.in,  get  bomb!*0,  and  user  has  to  either  provide  information 

c - about  element  X,  or  select  new  elements  (start  over). 

c - When  data  is  satisfactorily  known  about  X,  data  is  assigned  to  (now) 

c - known  element  “A”,  which  the  rest  of  the  program  uses.  Process 

c - repeats  for  unknown  X  — >  known  B 


CALL  getX (chA , ist A . HVA , DA2 , IA2 , TA2 , bomb) 

IF  (bomb.NE.O.)  THEN 
iok=l 

WRITE(*,*) 

HRITE(*,*) 'Oops!  Couldn”t  find  ’,chA,’  in  file  (ABX.in)’ 
WRITE(*,*) 'Make  a  new  element?  (l=yes) ’ 

HRITE(*,’(”  [”,Il.”]")’)iok 
READ(*,’(A)’)str 
IF  (str.NE.’  ’)  READ(str,*)iok 
IF  (iok.EQ.l)  THEN 


c - get  the  info  about  unknown  element  "X’’.  When  finished,  user 

c - has  option  of  appending  new  element  data  to  file  ABX.in,  so 

c - element  is  “known”  next  time  program  is  run. 

CALL  newX (chA , ist A , MWA , DA2 , IA2 , TA2) 


ELSE 

REWIND(UNIT=10) 

GOTO  300 
ENDIF 
ENDIF 

REWIND(UHIT=10) 

CALL  getX (chB , istB , MWB . DB2 . IB2 . TB2 , bomb) 

IF  (bomb.NE.O.)  THEN 
iok=l 

WRITEC*,*) 

WRITE(*,*)’0ops!  Couldn^'t  find  ’.chB,’ 
WRITE (*,*) ’Make  a  new  element?  (l=yes)’ 
WRITE(*,’(”  [”.Il.”]”)’)iok 
READ(*,’(A)’)str 
IF  (str.NE.’  ’)  READ(str.*)iok 
IF  (iok.EQ.l)  THEN 
CALL  newX (chB , istB . MWB . DB2 , IB2 , TB2) 
ELSE 
GOTO  300 
ENDIF 
ENDIF 

CLOSE (UNIT=10) 

IF  ((chA(2:2).EQ.’  ’) .AND. (chB(2:2) .EQ. ’  ’ 
name=chA (1:1) //chB (1:1) 

ELSEIF  (chA(2:2).Eq.’  ’)  THEN 
name=chA (1:1)// chB 
ELSEIF  (chB(2:2).EQ.’  ’)  THEN 
name=chA//chB(l : 1) 

ELSE 

name=chA//chB 

ENDIF 

OPEN  (UNIT=10.  FILE=’ABXY.in’.STATUS=’old’ 
CALL  getXY(name.RAB.IAB, TAB, bomb) 

IF  (bomb.NE.O.)  THEN 
iok=l 

WRITE(*,*) 

WRITE(*,*) ’Oops!  Couldn”t  find  ’.name,’ 

WRITE(*,*) ’Mcike  a  new  compound?  (l=yes)’ 

WRITE(*,’(”  [”,I1,”]  ”)’)iok 

READ(*,’(A)’)str 

IF  (str.NE.’  ’)  READ(str,*)iok 

IF  (iok.EQ.l)  THEN 


in  file  (ABX. in) ’ 


))  THEN 


) 


in  file  (ABXY.in)’ 
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c 

c 


c 

c 


c 

c 

c 


c 

c 

c 

c 


CALL  newXY(name,RAB.IAB.TAB) 

ELSE 
GOTO  300 
ENOIF 
ENDIF 

CL0SE(UNIT=10) 

WaiTEC* , ' (A) ’ ) ' +Using  ' //chA// ' . » //chB// ’.’//name 
WRITEC*,*) 

•get  ratios  of  element  A:B  (parts  of  A  vs.  parts  of  B) .  Converts 
•to  fractions  of  nuclei. 

WRITE(*,*) 'Stoichiometric  Proportions  of  ’//chA//’ , ’//chB//’?’ 
WRnE(*.’("  [".1PE11.3.’'.".1PE11.3.”]  ”)’)y0A,y0B 
READ(*.'(A)')str 

IF  (str.NE.'  ’)  READ(str,*)y0A,y0B 

WRITEC*. 192) '+Proport ions  of  '//chA//' , '//chB.yOA.yOB 
IF( (chA . EQ . ChB) . AND . (yOB . EQ . 0 . ) )  THEN 
y0A=0 . 
yOB=l . 

ENDIF 

WRITE(*.*) 

■Chamber  pressure  used  for  Saha  eqs.  Exit  pressure  used  for 


■ideal  rocket  calculations. 

WRITE(*,») 'Chamber,  exit  pressures  (atm)?' 

WRITE(*.'("  [",1PE11.3.",".1PE11.3."]  ”) ')pressl,press2 
READ(*,'(A)')str 

IF  (str.NE.'  ')  READ(str,*)pres5l,press2 
HRITE(*,192) '+pressl,press2  (atm) ' ,pressl,press2 
HRITE(*,*) 


-get  chamber  temperature  range  (independent  variable  iterated) 

-If  howmany*!,  program  operates  on  first  tempature  only,  regardless 
-of  upper  temperature. 

WRITE(*,*) 'Chamber  temperature  range:  tempi (K) ,temp2(K) , howmany?' 

WRITE(*,'("  C",1PE11.3,",",1PE11.3,",",I5, "]  ")') 
ttempl , temp2 , howmany 

READ(*,'(A)')str 

IF  (str.NE.'  ')  READ(str,*)templ,temp2, howmany 

WRITE(*, '(A,1P2(1PE11.3) ,15) ') 


A '  -t-tempi  (K) ,  t  emp2  (K) ,  howmany ' ,  tempi ,  t  emp2 ,  howmany 
-all  interesting  information  is  printed  to  file  if  running  only 
-one  temperature;  if  multiple  temperatures  are  run,  program  prints 
-output  data  to  sepaurate  files;  due  to  file  limitations  of  MS-DOS, 
-user  has  to  choose  between  molefraction  and  density  outputs.  Sigh... 


40 


IF  (howmany.GT. 1)  THEN 
WRITEC*,*) 

WRITE(*,*) 'Print  density  (1)  or  mole  fraction  (2)  ?’ 

WRITE(*.’("  [".II.”]  ")’)densmol 
READ(*.’(A)')str 

IF  (str.NE.’  ’)  READ(str,*)densmol 
WRITE(*.’("+densmol  [".II."]  ")')densmol 
ENDIF 
WRITEC*.*) 

c - On  personal  computers,  saha  equations  may  take  a  long  time  to  converge, 

c - especially  if  no  numeric  co-processor  is  used.  To  entertain  the  user, 

c - and  reassure  them  that  the  program  is  running  (instead  of  hanging), 

c - variables  seein.seemid,  and  maxout  count  the  number  of  times  the  inner, 

c - middle,  and  outer  nested  while  loops  proceed  before  pooping  out  to  the 

c - screen  the  loop  count  and  the  convergence.  If  the  computer  is  fast, 

c - printing  out  to  the  screen  every  loop  really  slows  the  program  down; 

c - but  if  the  computer  is  slow,  the  user  can't  tell  if  the  program  is 

c - working  or  not. 

c - Sugesstions:  fast  machines,  set  seein,seemid,maxout=100, 100,20 

c - slow  machines,  seein,seemid,maxouts20,20,20. 

c - The  program  is  "uone"  when  the  final  residual  (sum  of  all  the 

c - absolute  values  of  the  saha  equation  Ihs-rhs)  is  less  than  epsilon, 

c - where  epsilon=le-7.  When  the  outer  loop  residual  gets  down  around 

c - le-6,  it  starts  bouncing  around,  le-6  is  still  pretty  small,  so  if 

c - the  prograum  isn't  done  before  maxout  outer  loops,  the  program  stops. 

WRITE(*,*) 'View  loop  progress:  inner, middle, outer?' 

WRITE(*,'("  [".14. ",".14. ".",14,"]  ")') 
ftseein , seemid .maxout 
READ(*,'(A)')str 

IF  (str.NE.'  ')  READ(str,*) seein.seemid, maxout 
WRITE(*, 113) '+seein, seemid, maxout ' ,  seein, seemid, maxout 
WRITE(*.*) 

c - for  rocket  stuff 

WRITE(*,*) 'mass  flow  rate  (kg/s)?' 

WRITE(*,'("  [".1PE11.3. "]  ")')mdot 
READ(*, '(A)')str 
IF  (str.NE.'  ')  READ(str,*)mdot 
WRITE(*,111) '+mdot: ' ,mdot 
WRITE(*.*) 

REWIND(UNIT=1) 

c - write  info  to  file  AB.in  for  next  time  defaults 

WRITE (l,*)chA,chB 
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VJRITE(l,*)yOA,yOB 

WRITECl . *)pressl ,press2 

VRITECl ,  *)  tempi ,  temp2 ,  howmainy 

WRITE(l,*)deiismol 

VRITECl , *) seein . seemid .maxout 

VRITECl, ♦)mdot 

CL0SECUNIT=1) 

c - convert  atm’s  to  Pa’s 

pressl=pressl*atm 

pres82=press2*atm 

total=yOA+yOB 

yOA=yOA/total 

yOB=yOB/total 

c - open  appropriate  data  files  and  get  atomic  spectroscopic  data; 

c - load  spectroscopic  data  into  array  name  ''gu’’:  g  for  degeneracy, 

c - u  for  electronic  excitation  energy  level;  ionization  energies, 

c - in  cm“~l,  go  into  array  "ion” 

CALL  opunCchA.istA) 

CALL  getguCistA,chA.guA,ionA) 

CALL  cloze 

CALL  opunCchB.istB) 

CALL  getguCistB,chB,guB,ionB) 

CALL  cloze 

END  {subroutine  startup 

C*^^******^Hi^l^^^^^^^^^^^^^*t^^*^^**t******************************************** 
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A. 3  getX.for 

c — get  data  from  file  ABX.in  (unit  10)  for  unknown  element  X: 
c — using  element  neune  (ch),  gets  highest  element  state  to  consider,  i.e. 
c — highest  ionization  level  +  1  (ist),  gets  molecular  weight  (MW,  in 
c — g/mol),  gets  dissocation  (uD)  and  ionization  (ul)  energies  of  the 
c — diatomic, homonuclear  molecule,  in  eV,  and  gets  the  characteristic 
c — temperatures  of  rotation  and  vibration  for  the  neutral  and  singly 
c — ionized  diatomic  molecule  (in  K) .  If  cannot  find  ch  in  the  file, 
c — sets  bomb  to  non-zero,  which  flags  the  program  outside  this  routine. 
C********************************************************************** 
SUBROUTINE  getX( 

4ch , ist , MW , uD , ul , charT , bomb) 

CHARACTER+2  ch,readch 
INTEGER  ist 

REAL*8  MW,uD,uI,charT(4) ,bomb 
bomb=0 . 

READ (10,*) 

READdO ,  ♦ ,  ERR=100 ,  END=100)  readch ,  ist ,  MW,  uD ,  ul ,  charT 
DO  WHILE  (readch. NE.ch) 

READ ( 10 , * , ERR=100 , END=100) readch , ist , MW , uD , ul , charT 
ENDDO 
RETURN 
100  bomb=l. 

END 
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A. 4  getXY.for 

c — using  the  diatomic  molecule  name,  composed  of  the  compressed, 
c — concatenated  element  names,  in  alphabetical  order,  subroutine  looks 
c — in  file  ABXY.in  for  information  about  the  molecule.  Gets:  reaction 
c — (RAB)  and  ionization  (lAB)  energy  of  diatomic,  heteronuclear 
c — molecule,  and  the  characteristic  temperatures  for  rotation  and 
c — vibration  for  the  neutral  and  singly  ionized  molecule  (both  at  ground 
c — electronic  state).  If  cannot  find  info,  bomb!=l.  If  molecule  is 
c — pseudo-heteronuclear,  e.g.  it  looks  for  diatomic  nitrogen,  even 
c — though  N2  is  homonuclear,  homonuclear  info  is  returned. 

SUBROUTINE  getXY( 
tname , RAB , I AB , TAB , bomb) 

CHARACTER*4  name,readname 
REALMS  RAB , I AB , TAB (4) , bomb 
bomb=0 . 

READ (10,*) 

READ (10 , ♦ , ERR=100 , END=100) readname , RAB , lAB , TAB 
DO  WHILE  (readname. NE. name) 

READ ( 10 , * , ERR=100 , END= 100) readname , RAB , I AB , TAB 
ENDDO 
RETURN 
100  bomb=l. 

END 

Q********^l*^^***^*^^t*^t^^***:i^HHH^1^***************^^************************* 


44 


A.  5  newX.for 


c — Prompts  user  for  information  about  new  elements.  Asks  for 
c — information  that  would  otherwise  be  found  in  file  ABX.in.  After 
c — completion,  gives  user  opportunity  to  append  new  element  info  to  file 
c — ABX.in,  so  next  time  the  information  will  be  found.  Gets  highest 
c — ionization  state  (ist),  molecular  weight  (MW),  dissocation  (uD)  and 
c — ionization  (ul)  energies  of  diatomic  homonuclear  molecule,  in  eV,  and 
c — characteristic  temperatures  of  rotation  and  vibration,  at  the  ground 
c-electronic  state,  of  the  neutral  and  singly  ionized  diatomic  molecule. 

SUBROUTINE  newX( 
tch, ist, MW, uD,uI, chart) 

CHARACTER  ch*2,str*30 
INTEGER  ist 

REAL*8  MW,uD,uI,charT(4) , iok 
10  WRITEC*,*) ’Let’ 's  mzdce  a  new  element...’ 
ist=0 
MW=0. 
uD=0. 
ul=0. 

charT(l)=l. 
charT(2)=l. 
charT(3)=l . 
charT(4)=l. 

WRITEC*,*) 

WRITE(*,*) ’Highest  atomic  state  (1+highest  ienization  state)?  ’ 

WRITE(*,’(”  [”,I1,”]  ”)’)ist 

READ(*,’(A)’)str 

IF  (str.NE.’  ’)  READ(str,*)ist 

WRITE(*,’(”+ist:  [”  ,11,  ”]  ”) ’)ist 

WRITEC*,*) 

WRITEC*,*) ’Atomic  mass  (amu)?  ’ 

WRITE(*,’(”  [”,1PE11.3,”]  ”)’)MW 
READ(*,’(A)’)str 
IF  (str.NE.’  ’)  READ(str,*)MW 
WRITE(*,’(”+MW:  [”  ,1PE11.3,  ”]  ”) ’)MW 
WRITE(*,*) 

WRITE(*,*) ’Dissociation  energy  of  diatomic,  homonuclear  ’// 
ft’molecule  (eV)  ?’ 

WRITE(*,’(”  [”,1PE11.3,  ”]  ”)’)uD 

READ(*,’(A)’)str 

IF  (str.NE.’  ’)  READ(str,*)uD 


WRITE(*,’("+uD:  [' MPE11.3,  ’ »]  ”)  OuD 
WRITEC*,*) 

WRITEC*,*) 'Ionization  energy  of  diatomic,  homonuclear  ’// 

R 'molecule  (eV)  ?' 

WRITE(*.'("  [".1PE11.3."]  ")')ul 
READ(*,'(A)')str 
IF  (str.NE.'  ')  READ(str,*)uI 
WRITE(*.'("+uI:  C".lPE11.3."]")')uI 
WRITEC*,*) 

HRITEC*,*) 'Chajracteristic  temperatures,  at  ground  electronic’// 
t' state  (K)  ?' 

WRITEC*,*) 'Neutral:  Rotation,  Vibration;  Singly  Ionized:  ’// 
^'Rotation,  Vibration' 

WRITE(*,'("  [",lP4(E11.3),"]")')charT 

READ(*, ' (A) ')str 

IF  (str.NE.'  ')  READ(str,*)charT 

WRITEC* , ' (  "  +charT :  [  "  , 1P4(E1 1 . 3) .  " ]  " ) ' ) charT 

WRITEC*,*) 

WRITE(*,*) 'Summary: ' 

WRITEC*,*) 'Highest  atomic  state  (1+highest  ionization  state)  ' 
WRITE(*,'("  [".II,"]  ")')ist 
WRITEC*,*) 'Atomic  mass  (amu)  ' 

WRITE(*.'("  [",1PE11.3,"]")')MW 

HRITEC*,*) 'Dissociation  energy  of  diatomic,  homonuclear  ’// 
ft'molecule  (eV)  ' 

WRITEC*, '  ("  [' '  .1PE11.3,  "]  ")  ')uD 

WRITE(*,*) 'Ionization  energy  of  diatomic,  homonuclear  '// 
t'molecule  (eV)  ' 

WRITEC*, '("  [",lPE11.3."]")')uI 

WRITEC*,*) 'Characteristic  temperatures,  at  ground  electronic’// 
ft 'state  (K)  ' 

WRITEC*,*) 'Neutral:  Rotation,  Vibration;  Singly  Ionized:  ’// 
ft 'Rotation,  Vibration' 

WRITEC* ,  ’  (  "  ["  .  1P4(E11 .3) ,  "]  " ) ' )chaiT 
iok=l 

WRITEC*,*) 'Is  everything  ok?  (l=ok,other=redo) ' 

WRITEC*, '("  [".Il."]")')iok 

READC*, '(A) ')str 

IF  (str.NE.'  ')  READ(str,*)iok 

IF  (iok.NE.l)  GOTO  10 

iok^O 

WRITEC*,*) 'Good.  Append  data  to  file  (ABX.in)?  (0=no , l=yes) ’ 


WRITE(*,’('»  [".11,"] 'OOiok 
READ(*.'(A)’)str 
IF  (str.NE.'  ')  READ(str.*)iok 
IF  (iok.EQ.l)  THEN 

WRITEC*,*) ’Appending  data  to  file  (ABX.in)...’ 
WRITE(10.’(A,I2.1P7(E11.3))’)  ch. ist.MW.uD.uI.charT 
ELSE 

WRITEC*,*) ’Will  not  append  data  to  file...’ 

ENDIF 

END  'subroutine  newX 
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A. 6  newXY.for 


c — Asks  user  for  information  about  unknown  compound.  Gets  reaction 
c — (RAB)  and  ionization  (lAB)  energy  of  heteronuclear  diatomic  molecule, 
c — in  eV.  Also  gets  characteristic  temperatures  of  rotation  and 
c — vibration,  at  the  ground  electronic  state,  of  the  diatomic  neutral 
c — and  singly  ionized  molecule.  User  can  append  data  to  end  of  input 
c — file  ABXY.in  to  save  typing  data  in  again  next  time. 

C********************************************************************** 
SUBROUTINE  neuXYC 
kname,RAB, IAB,TAB) 

CHARACTER  name*4,str*30 
REAL*8  RAB,IAB,TAB(4) 

WRITE(*,*) 

10  WRITE(*,*) 'Lef’s  madse  a  new  compound...’ 

RAB=0. 

IAB=0. 

TAB(1)=1. 

TAB(2)=1. 

TAB(3)=1. 

TAB(4)=1. 

WRITEC*,*) 

WRITEC*,*) 'Reaction  energy  of  diatomic,  heteronuclear  ’// 
t ’molecule  (eV)  ?’ 

WRITE(*,'("  C’MPE11.3,"]'0')RAB 

READ(*,’(A)’)str 

IF  (str.NE.’  ’)  READ (str,*) RAB 

WRITEC*,  ’(”+RAB:  [”  ,1PE11.3,  ”]  ”) ’)RAB 

WRITEC*,*) 

WRITEC*,*) ’Ionization  energy  of  diatomic,  heteronuclear  ’// 
ft ’molecule  (eV)  ?’ 

WRITEC*, ’(”  [”,1PE11.3,”]  ”)’)IAB 

READ(*,’(A)’)str 

IF  (str.NE.’  ’)  READ(str,*)IAB 

WRITEC*, ’(”+IAB;  C”,1PE11.3,”]  ”)’)IAB 

WRITEC*,*) 

WRITEC*,*) 'Characteristic  temperatures,  at  ground  electronic’// 
ft 'state  (K)  ?’ 

WRITEC*,*) 'Neutral:  Rotation,  Vibration;  Singly  Ionized:  ’// 
ft 'Rotation,  Vibration’ 

WRITEC*, ’(”  [”,1P4(E11.3),  ”]  ”)’)TAB 

READ(*,’(A)’)str 

IF  (str.NE.’  ’)  READ(str,*)TAB 
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WRITE(»,'(”+TAB:  [' MP4(E11 .3) .  * ’  O ’)TAB 
iok=l 

WRITEC*,*) 

WRITEC*,*) 'Reaction  energy  of  diatomic,  heteronuclear  '// 
ft 'molecule  (eV)  ' 

HRITE(*.'("  [".1PE11.3,"]  ")’)RAB 

WRITEC*,*) 'Ionization  energy  of  diatomic,  heteronuclear  ’// 
ft 'molecule  (eV)  ' 

WRITEC*, '("  [",1PE11.3."]  ")')IAB 

WRITEC*,*) 'Characteristic  temperatures,  at  ground  electronic’// 
ft 'state  CK)  ' 

WRITEC*,*) 'Neutral:  Rotation,  Vibration;  Singly  Ionized;  '// 
ft’Rotation,  Vibration’ 

WRITEC* , ’ C ’ '  [ ”  , 1P4CE11 . 3) .  ”]  ”  ) ’ )TAB 
WRITEC*,*) 'Is  everything  ok?  Cl=ok,other=redo) ’ 

WRITEC*, ’C”  [”.Il.”]”)’)iok 

READC*,’CA)’)str 

IF  Cstr.NE.’  ’)  READCstr.*)iok 

IF  Ciok.NE.l)  GOTO  10 

iok=0 

WRITEC*,*) 'Good.  Append  data  to  file  CABXY.in)?  C0=no, l=yes) ’ 

WRITEC*, ’C”  C".Il,”]")’)iok 

READC*.’CA)’)str 

IF  Cstr.NE.’  ')  READCstr,*)iok 

IF  Ciok.EQ.l)  THEN 

WRITEC*,*) 'Appending  data  to  file  CABXY.in)...’ 

WRITECIO, ’CA,lP6CE11.3))’)name,RAB,IAB.TAB 
ELSE 

WRITEC*,*) 'Will  not  append  data  to  file...’ 

ENDIF 

END 

C********************************************************************** 
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A. 7  opun.for 

c — open  the  necessary  files  to  get  the  atomic  spectra  data.  The  naming 
c — convention  is  as  follows:  the  input  file  names  start  with  ch  (removed 
c — of  spaces),  have  an  underscore  and  then  an  number,  0-4, 

c — which  corresponds  to  the  charge  on  the  atom.  Charge  -<■4  corresponds 
c — to  ist=5,  e.g.  C_V,  carbon  ionized  four  times. 

SUBROUTINE  opun( 
tch,i8t) 

CHARACTER  ch*2,filenum*l, filename* 12 
INTEGER  i 
DO  10  i=0,ist-l 
WRITECfilenum, *(Il)»)i 
0PEN(UNIT=10+i. 

k  FILE=filename(ch, ’_*//f ilenum//’ .in'), 
k  STATUS='0LD') 

10  CONTINUE 

END  ! subroutine  opun 

Q^^***^^l^^^^^l:^^^^*^t^*^^*l|^^t*^*^^:^t^^i^i^^i^:^i*^^i^i|^^i*^i**^i********************************* 


A. 8  cloze.for 

c — close  all  the  opened  files 

Q*^^*^^t^^^^^^*^^*iilti^i^^***:t^**^H^**^^^^:tli^^^^^^:^^l|^^^^^t********************************** 

SUBROUTINE  cloze 
CL0SE(UNIT=10) 

CL0SE(UNIT=11) 

CL0SE(UNIT=12) 

CL0SE(UNIT=13) 

CL0SE(UNIT=14) 

EIND  !  subroutine  opun 

Q**:tf:ilfJ^^f*ii**ilf*tt********it*****:li***i^**:************************************ 
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A.9  getgu.for 

c — necessary  files  are  already  opened,  so  reads  from  files  unitslO 
c-'through  unit*10+ist-l  to  get  the  ionization  energy  of  the  atom,  in 
c— cm‘-l,  and  the  (degeneracy,  energy  level)  pairs  (energy  levels  also 
c — in  cm“-l).  Ionization  energies  stored  in  array  “ion'\  and 
c — degeneracy,  level  pairs  stored  in  array  “gu”.  For  example, 
c — guA(5,200,2)  stores  the  200x2  (g,u)  pairs  for  the  5  species  of 
c — element  A. 

C*^i******t*****************************ilf****************************^t** 

SUBROUTINE  getgu( 

Aist,ch,gu,ion) 

INTEGER  ist 

CHARACTER  ch*2,f ilenum’fl ,f ilename’^12 
REALMS  gu(5.200,2) ,ion(5) 

INTEGER  i.j.k 
IMPLICIT  REAL*8  (a-h,l-z) 

DO  100  i=1.5 
DO  100  j=1.200 
DO  100  k=1.2 

100  gu(i,j.k)=0. 

DO  200  i=l,ist 
WRITE(filenum. '(Il)‘)i-1 
WRITE(*,*) ’Getting  data  from  file  ’// 

A  filenajne(ch, ’_’//filenum//’-in*) 
idisk®9+i 

READ(idisk,«)  ion(i) 
j=0 

READ(idisk,*)  degen, englev 
DO  WHILE  (degen. GT.O.) 

j=j+l 

gu(i,j,l)=degen 
gu(i,j,2)*englev 
READ(idisk,*)  degen, englev 
ENDDO 

200  CONTINUE 

END  {subroutine  getgu 

C*******************<t‘*********************4‘**************************** 
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A. 10  saha.for 


c — The  meat  of  the  program.  Subroutine  has  3  nested  while  loops,  one 
c — of  which  is  hidden  in  SUBROUTINE  balance.  Outer  loop  updates 
c — electron  and  heavy  fractional  densities  (ye,  yh) ,  and  beta 
c — (beta=ye/(yh+ye)) .  Middle  loop  holds  ye, yh, beta  constant,  but 
c — updates  heteronuclear  reaction  products  AB,  AB+.  Inner  loops  update 
c — atomic  and  homonuclear  diatomic  species.  (See  SUBROUTINE  balance  for 
c — inner  loops.) 

SUBROUTINE  saha( 

ftpressl , temp, istA, istB , seein , seemid .maxout , 

&MWA . DA2 , I A2 . TA2 . ionA , guA . 

&MWB , DB2 , IB2 , TB2 , ionB , guB , 
ftZOA , ZA , ubar A , uubA , rhsA , yOA , yA , 

&Z0B , ZB , ubaxB , uubB , rhsB , yOB , yB , 
feenthA , entrA , CpA , 
ftenthB , entrB , CpB , 

ftRAB , I AB , TAB , ZOAB , ZAB , ubar AB , uubAB , rhs AB , yAB , 
tenthAB , entr AB , CpAB , enthe , Cpe , entre , ye , yh , beta) 

IMPLICIT  REALMS  (a-h,l-z) 

INTEGER  ist A , istB , see in , seemid , maxout 
REALMS  pressl.temp, 

4MWA , DA2 , IA2 . TA2 (4) . ionA (5) . guA (5 , 200,2), 
ftMWB,DB2,IB2,TB2(4) , ionB(5) ,guB(5,200,2) , 
ftZ0A(7,3) ,ZA(7,3) ,ubarA(7) ,uubA(7) ,rhsA(7) ,yOA,yA(7) , 
4Z0B(7,3),ZB(7,3),ubarB(7),uubB(7),rhsB(7),y0B,yB(7), 

4enthA(7) ,entrA(7) ,CpA(7) , 

4enthB(7) ,entrB(7) ,CpB(7) , 

4RAB,IAB,TAB(4), 

4Z0AB(2,3) ,ZAB(2,3) ,ubarAB(2) ,uubAB(2) ,rhsAB(2) ,yAB(2) , 

4enthAB(7) ,entrAB<7) ,CpAB(7) , 

4ye,yh,beta 

C  local  variables  for  subroutine  saha 

INTEGER  i,imid,iout,ityA,ityB 
REAL*8 

4epsilon , chi , ymin , 

4yeA,yhA,resA, 

4yeB,yhB,resB, 

4res AB , wyh , residual 

DATA  epsilon, chi, ymin  /l.E-7,0.5,l.E-26/ 

c - get  the  partition  functions  for  species  of  A,B,  and  AB  for  the 

c - given  temperature  and  pressure. 


53 


CALL  Z1 (temp . istA , guA . TA2 . lonA , ZA , ubar A , uubA) 

CALL  Z1 (temp , istB . guB . TB2 , ionB , ZB . ubarB , uubB ) 

CALL  Z2 (temp . ZA , ubar A , uubA , ZB , ubarB , uubb , TAB , ZAB , ubar AB , uub AB ) 

c - get  the  right-hand-sides  (rhs)  of  the  saha  equations  (R_i  in  the 

c - report),  which  are  specified  for  a  given  temperature  and  pressure. 

CALL  rhsl (temp , pressl , istA , ionA , ZA . DA2 , IA2 , MWA , rhsA) 

CALL  rhs 1 (temp , pressl . istB , ionB , ZB , DB2 , IB2 , MWB , rhsB ) 

CALL  rhs2 (temp , pressl . ZA . ZB . ZAB , RAB , I AB . HVA , HUB . rhs AB) 

c - get  thermodynamic  properties  for  the  individual  species;  depend 

c - only  on  temp,  press,  pajrtition  functions,  and  ubar,  uub;  can  get 

c - whole  gas  thermodynamic  properties  later  by  tzdcing  the  weighted 

c - average  thermodynamic  properites,  where  the  weights  are  the  mole 

c - fractions. 

CALL  therml (teoq) , pressl , DA2 , IA2 , ionA , ZA , ZO A , ubar A , uubA , 

AenthA , CpA , entrA) 

CALL  therml (temp , press 1 , DB2 , IB2 , ionB , ZB , ZOB , ubarB , uubB , 
henthB . CpB . entrB) 

CALL  tberm2 (temp , pressl , DA2 , 0B2 , RAB , I AB , ZAB , ZOAB . 

AubarAB , uubAB , enthAB , CpAB , entrAB) 

CALL  therme (temp , pressl , enthe , Cpe . entre ) 

DO  200  i=1.7 
yA(i)=0. 
yB(i)»0, 

200  CONTINUE 

DO  300  i=1.2 
yAB(i)=0. 

300  CONTINUE 

c - guess  (“ges”)  the  initial  values  of  the  atomic  and  homonuclear 

c - atomic  species  fractions  by  assuming  that  (1)  elements  A,B  can 

c - mix,  but  do  not  react  chemically,  (2)  a  given  element  has  two 

c - dominant  species  present  at  any  one  time,  all  others  being 

c - negligibly  (sp?)  small.  For  a  single  element,  or  for  mixtures, 

c - this  method  will  give  the  right  answer  to  about  57.,  and  zippy 

c - quick  too!  But  chemical  reactions  between  A  and  B  require  less 

c - elegant  methods . . . 

CALL  gesl(rhsA,yOA,yA) 

CALL  ge8l(rhsB,y0B,yB) 

c - itery  finds  the  dominant  species  in  an  element,  and  sets  the 

c - variably  ity  to  that  species  position  in  the  element  array.  In 

c - subsequent  calculations,  operations  are  performed  using  that 

c - species  as  the  “row  pivot”.  This  reduces  numerical  error, 

c - which  means  the  set  of  equations  converges  more  quickly  to  a  good 
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c - answer. 

CALL  itery(yA,ityA) 

CALL  itery(yB,ityB) 

c - functions  that  get  the  sum  of  the  atomic  fractions  for  the 

c - “heavy”  fraction  (yh) ,  and  the  weighted  average  atomic 

c - fractions  for  the  electron  fraction  (ye),  where  the  weights  are 

c - the  species  charges.  At  this  point  in  the  program,  this  is  just 

c - our  first  guess  at  ye,yh,and  beta. 


CALL  getye(yA,yeA) 

CALL  getye(yB,yeB) 

CALL  getyh(yA,yhA) 

CALL  getyh(yB,yhB) 
ye=yeA+yeB+yAB(2) 
y h=yhA+yhB+yAB ( 1 ) +yAB ( 2 ) 
beta=ye/(yh+ye) 


c - calculate  the  maximum  possible  fraction  of  species  AB,  as  if  all 

c - of  one  species  or  another  combined  to  form  the  compound;  maximum 

c - is  modified  by  the  rhsAB,  which  sort  of  indicates  how  much  the 

c - species  AB  dissociates  as  the  temperature  increases. 


ymaxAB=0 . 5-DSqRT(0 . 25-yOA*yOB/ ( 1 . +rh8 AB ( 1 ) ) ) 

iout=0 

residual*!. 


c - outer  while  loop:  finished  when  either  (1)  residuaKepsilon,  which 

c - means  the  y’s  we  have  satisfies  the  system  of  equations  nicely, 

c - or  (2)  the  loop  has  tried  maxout  times  to  converge,  and  hasn’t 

c - done  it  yet,  and  so  probably  won’t,  maxout  is  definable  by  the 

c - user  in  the  input  section:  maxout=20  works  well,  and  reasonably 

c - balemces  speed  vs.  efficiency.  This  loop  balamces  charge. 


DO  WHILE  ((residual. GT. epsilon) .AND. (iout.LT. maxout)) 
yeO=ye 
yhO=yh 
imid=0 

c - reinitialize  the  upper  and  lower  bounds  for  yAB 

yhiAB=ymaxAB 
yloAB=0. 
wyh=-i . 


c - middle  while  loop:  is  done  when  the  weighted  sum  of  the  heavy 

c - atomic  fractions  (wyh)  is  within  epsilon  of  1,  where  the 

c - weights  are  the  number  of  nuclei  in  the  molecule  (2  for 

c - diatomic,  1  for  atomic).  This  is  essentially  the  conservation 

c - of  mass  equation.  Within  this  loop,  charge  is  held  temporarily 

c - constant. 
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DO  VHILE  (DABS(vyh-l.).GT. epsilon) 


yABO»yAB(l) 

c - - - inner  loopCs),  one  each  for  the  elements  A  and  B:  is  done  when 

c - —the  weighted  sum  of  the  atomic  fractions  (sum  of  element 

c - —nuclei)  is  within  epsilon  of  the  total  atomic  fraction,  yO; 

c - - — This  is  essentially  conservation  of  mass  for  each  element. 

CALL  balance (ye . yh . beta, epsilon, chi , ymin , yAB , seein , 
k  yOA,yA,ityA,rhsA) 

CALL  balance (ye , yh , beta, epsilon , chi , ymin , yAB .seein . 
k  yOB,yB,ityB,rhsB) 

c - done  with  inner  loops,  so  now  update  yAB,  yAB+, 

c - -under-relaxing  using  bisection  method.  Converges  quickly  for 

c - monotonic  functions. 


yAB ( 1 ) =yA ( 1 ) ♦yB ( 1 ) / (yh+ye) /rhsAB ( 1 ) 
yABl=yAB(l) 

IF  (yABl.GT.yABO)  THEN 
yloAB=yAB0 

yAB (l)=(l.-chi) *yABO+chi*yhiAB 
ELSEIF  (yABl.LT.yABO)  THEN 
yhiAB=yABO 

yAB(l)=(l , -chi)*yABO+chi*yloAB 
ENDIF 

IF  (ye.EQ.O.)  THEN 
yAB(2)=0 

ELSEIF  (ye.GT.O.)  THEN 
yAB (2) =yAB ( 1 ) *rhs AB (2) /beta 
ENDIF 

c - get  new  weighted  sums  of  heavies  for  middle  loop  check 

CALL  getwyh(yA,yAB,wyhA) 

CALL  getwyh(yB,yAB,wyhB) 

wyh=wyhA+wyhB 

imid=imid+l 


922  FORMAT  (’+’ ,22x,lP2(E11.3)) 

c - the  progress-o-meter 

IF  (mod(imid,seemid).EQ.O)  THEN 
WRITE(* . 922) imid , wyh-1 . 

ENDIF 


ENDDO 

c - get  new  electron,  heavy  fractions  for  the  outer  loop  check 

CALL  getye(yA,yeA) 

CALL  getye(yB,yeB) 

CALL  getyh(yA,yhA) 


56 


c 

c 

c 


c 


944 

c— 


CALL  getyh(yB,yhB) 

ye=yeA+yeB+yAB(2) 

yh=yh A+y hB+y AB (l)+yAB(2) 

yel=ye 

yhl=yh 

•update  electron,  heavy  density  fractions  by  under-relaxing, 

•using  bisection  method.  Again,  best  for  monotonic  functions 

•(such  as  ve  have). 

ye=(l . -chi)*yeO+chi*yel 

yh=(l . -chi)*yhO+chi*yhl 

beta=ye/(yh•^ye) 

■get  partial  and  total  residuals  for  outer  loop  check 
CALL  resl (yOA , yA , rhsA , yAB , ye , yh , bet a , res A) 

CALL  resl (yOB , yB , rhsB , yAB , ye , yh , beta , resB) 

CALL  re82 (yA , yB , yAB , rhsAB , ye , yh , beta , res AB) 

residual=resA-^resB■•■resAB 

iout=iout+l 

FORMAT  (H',44x,1P2(E11.3)) 
are  we  almost  there? 

WRITE(* , 944) iout , residual 


ENDDO 

END  ! subroutine  saha 
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A.ll  Zl.for 

c — Make  the  partition  function  for  the  element,  given  the  atomic  spectra 
c — data  in  arrays  gu  and  ion,  and  the  temperature.  Also  calculates  the 
c — average  electronic  excitation  energy  level,  ubar,  and  the  average 
c — squared  energy  excitation  energy  level,  uub  (short  for  (u*u)bar). 
C"These  are  useful  later  in  calculateing  thermodynamic  properties,  such 
c— as  heat  capacity.  Array  ZA(7,3),  holds,  for  example,  partition 
c-f unction  data  for  the  7  states  (atomic  neutral,  atomic  +l,...+4, 
c — diatomic  neutral,  diatomic  +1)  of  element  A,  for  electronic, 
c — rotational,  and  vibrational  excitation  (Zrot  and  Zvib  =1  for  atomic 
c — species).  Electronic  partition  functions  for  diatomic  species  are 
c — approximated  by  appropriately  averaged  atomic  partition  functions, 
c — For  exsusple,  Zel,A2-*-  Xapprox  (ZelA)*(ZelA+) .  Remember,  partition 
c — functions  multiply,  not  add. 

SUBROUTINE  Zl( 

ftt emp,i8t,gu, chart, ion, Z, ubar, uub) 

INTEGER  ist 

REAL*8  temp,gu(5,200,2),charT(4),ion(5),Z(7,3),ubar(7) ,uub(7) 
INTEGER  i,j 

IMPLICIT  REAL*8  (a-h,l-z) 
tempinv=l . 4387/temp 
DO  100  i=l,7 
DO  100  j=l,3 
100  Z(i,j)=0 

DO  200  i=l,ist 
summ=0 . 
sumu=0 . 
sumuu=0 . 

j=l 

DO  WHILE  (gu(i,j,l).GT.0.) 
term  =  0. 

IF((gu(i, j,2)*tempinv.LT.70.) 
k  .AND. (gu(i, j ,2) .LT.ion(i)-temp/1.4387)) 
k  term=gu(i,j,l)*DEXP(-gu(i,j,2)*tempinv) 

summ  =^8umm  't’term 
sumu  =sumu  +term*gu(i, j ,2) 

8umuu=sumuu+term*gu(i, j ,2)*gu(i, j ,2) 

ENDDO 

Z(i,l)=summ 

Z(i,2)=l. 
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Z(i.3)=l. 
ubar ( i ) =sumu/summ 
uub  ( i) =suffluu/siunm 
200  CONTINUE 

Z(6.1)=Z(1,1)*Z(1,1) 

Z  (  6 , 2)  =temp/ chairT  (1) 

Z(6,3)=temp/charT(2) 
ubar (6 ) =ubar ( 1 ) +ubar ( 1 ) 

uub(6)=uub(l)-ubar(l)*ubar(l)+uub(l)-ubar(l)*ubar(l) 

Z(7.1)=Z(1.1.)*Z(2.1) 

Z(7,2)=temp/charT(3) 

Z(7 , 3) =temp/ charT (4) 
ubar(7)=ubar(l)+ubar(2) 

uub(7)=uub(l)-ubar(l)*ubar(l)+uub(2)-ubar(2)*ubar(2) 
END  {subroutine  Z 
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A. 12  Z2.for 


c — Gets  partition  functions  for  diatomic  species  (hence  the  2  in  Z2) . 
c — Approximates  the  electronic  partition  function  by  appropriately 
c — average  atomic  partition  functions.  For  example,  Zel,AB+  is  an 
c—average  between  Zel.A,  Zel,B+,  Zel,A+,  and  Zel.B.  Although  this  is 
c — not  rigorously  correct,  it  gets  in  the  ballpark;  I  don’t  worry  too 
c — much  about  diatomic  species  since  (1)  this  code  is  meant  to  be  run  at 
c — high  temperatures,  where  there  aren’t  many  diatomic  species,  and  (2) 
c — Selph’s  Isp  code  does  a  great  (and  better,  too)  job  at  lower 
c — temperatures.  Also  gets  averaged  ubar  and  (u*u)bar  for  diatomic 
c — species.  Note  that  of  ZA=ZB,  ZAB=ZA2=ZB2,  as  should  be.  Thus, 
c — entering  in  two  identical  species  does  not  confuse  the  code. 
C********************************************************************** 
SUBROUTINE  Z2( 

fttemp . ZA , ubarA , uubA , ZB , ubarB . uubb , TAB , ZAB , ubarAB , uubAB) 

REAL’t'8  temp, 

kZA (7.3). UbarA (7) , uubA (7) , 
tZB(7,3) .ubarB(7) .uubB(7) , 
mB(4)  .ZAB(2,3)  .ubarAB(2)  .uubAB(2) 

IMPLICIT  REAL*8  (a-h,l-z) 

ZAB(1.1)=ZA(1.1)*ZB(1.1) 

ZAB(1.2)=temp/TAB(l) 

ZAB(l,3)=temp/TAB(2) 
ubarAB ( 1 ) =ubar A ( 1 ) +ubarB ( 1 ) 

uubAB ( 1 ) *uubA ( 1 ) -ubarA (1) *ubar A ( l)+uubB ( 1 ) -ubarB ( 1 ) *ubarB ( 1 ) 
ZAB(2.1)=0.5*(ZB(2.1)*ZA(1,1)+ZA(2,1)*ZB(1,1)) 

ZAB(2,2)=temp/TAB(3) 

ZAB(2.3)=temp/TAB(4) 

ubarAB (2) »0 . 5* (ubarA ( 1 ) +ubar A (2 ) +ubarB ( 1 ) +ubarB (2 ) ) 
uubAB (2) =0. 5* ( 

AuubA ( 1 ) -ubarA ( 1 ) ♦ubarA ( 1 ) +uubB ( 1 ) -ubarB ( 1 ) ♦ubarB ( 1 ) + 
ftuubA(2)-ubarA(2)^ubarA(2)+uubB(2)-ubarB(2)^ubarB(2)) 

END  {subroutine  Z2 

Q*i^*******^m*4f*************i^******************************************* 
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A. 13  therml.for 


c — gets  the  thermodynamic  functions  (enthalpy,  heat  capacity,  and 
c — entropy  difference  from  Si'P)  for  cin  element’s  species,  and  sticks 
c — them  in  ajrrays  enth(7),  Cp(7),  and  entr(7).  Uses  the  element’s 
c — partion  functions  at  STP  (ZO)  and  the  current  T,P  (Z) .  See  Theory 
c — in  paper  for  derivation  of  thermodynamic  functions. 

SUBROUTINE  thermlC 

fttemp , press 1 , uD , ul , ion , Z , ZO , ubar , uub , enth , Cp , entr ) 

REAL*8  temp,pressl,uD,ion(5) ,Z(7,3) ,Z0(7,3) , ubar (7) , uub (7) , 
&enth(7) ,CpC7) ,entr(7) 

INTEGER  i 

IMPLICIT  REAL*8  (a-h,l-z) 

DATA  gasR, cminv , atm/8 . 31439 , 8 . 06549E3 , 0 . 101325e6/ 

DATA  tmpC  /2.9815E2/ 
tempinv=l . 4387/temp 
dision=0.5*uD*cminv 
DO  100  i=l,5 

IF  (Z(i,l).EQ.O)  THEN 
enth(i)=0. 

Cp(i)=0. 

entr(i)=0. 

ELSE 

enth(i)=gasR*temp*(2.5+(dision+ubar(i))*tempinv) 

Cp ( i ) =gasR* (2.5+ (uub (i)-ubar(i) *ubar ( i ) ) *temp inv* t emp in v) 
entr(i)=gasR*(2.5*DL0G(temp/tmp0)+ubar(i)*tempinv+ 

&  DLOG (Z (i , 1 ) /ZO ( i , 1 ) ) -DLOG (pressl/atm) ) 
dision=dision+ion(i) 

ENDIF 

100  CONTINUE 

enth(6)=gasR*temp*(4.5+ubar(6)*tempinv) 

Cp(6)=gasR*(4.5+(uub(6)-ubar(6)*ubar(6))*tempinv*tempinv) 

entr(6)=gasR*(4.5*DL0G(temp/tmp0)+ubar(6)*tempinv+ 

&DL0G(Z(6,l)/Z0(6,l))-DL0G(pressl/atm)) 

enth(7)=gasR*temp*(4.5+(uI*cminv+ubar(7))*tempinv) 

Cp(7)=gasR*(4.5+(uub(7)-ubar(7)*ubar(7))*tempinv*tempinv) 

entr(7)=gasR*(4.5*DL0G(temp/tmp0)+ubar(7)*tempinv+ 

&DL0G (Z (7 , 1 ) /ZO (7 , 1 ) ) -DLOG (press 1/atm) ) 

END 

C +*  + *  +  +  *  +  +  +  +  +  ***  ^  + 1)1  +  4c  :4c  4e 1)^ ;tc  :|E ^  4i 4c  4c ^ :)i itt  1)1  itc  itcitc  * t  itc  It: *  t  :tc 


A. 14  therin2.for 


c — calculates  the  thermodynamic  functions  (h,c_p,s-sO)  for  the 
c — heteronuclear  diatomic  species.  Sticks  them  in  arrays  enthAB(2), 
c— CpAB(2),  and  entrAB(2). 

SUBROUTINE  therm2( 

ktemp , pressl . DA2 , DB2 . RAB , lAB . ZAB , ZOAB , 
kubar AB , uubAB , enthAB , CpAB , entr AB ) 

REAL*8  temp.pressl.DA2,DB2,RAB.ZAB(2.3) ,Z0AB(2,3) , 
ftubarAB(2) ,uubAB(2) ,enthAB(2) ,CpAB(2) , entrAB(2) 

IMPLICIT  REALMS  (a-h.l-z) 

DATA  gasR, cminv, atm/8 . 31439 ,8 .06549E3 ,0 . 101325e6/ 

DATA  tmpO  /2.9815E2/ 

tempinv=l .4387/temp 

dision=(0 .5*(DA2+DB2)-RAB)*cminv 

enthAB ( 1 ) =gasR*temp* (4 . 5+ (dis ion+ubar AB ( 1 ) ) *temp inv) 
CpAB(l)=gasR*(4.5+(uubAB(l)-ubarAB(l)*ubarAB(l) )*tempinv*tempinv) 
entrAB (1) =gasR* (4 . 5*DL0G (t emp/ tmpO)+ubar AB ( 1 ) *t empinv+ 
ftDLOG (ZAB (1,1) /ZOAB (1.1)) -DLOG (press 1/atm) ) 
dision=(0 . 5*(DA2+DB2)-RAB+IAB)*cminv 
enthAB (2) =gasR*temp* (4 . 5+ (dis ion+ubar AB (2) ) *t emp inv) 

CpAB(2)=gasR* (4 . 5+(uubAB(2)-ubarAB(2) *ubarAB(2) ) *tempinv*tempinv) 
entrAB(2)=gasR*(4.5*DL0G(temp/tmp0)+ubarAB(2)*tempinv+ 

4DL0G (ZAB (2.1) /ZOAB (2,1)) -DLOG (press 1/atm) ) 

END 

C****^!**************************  ****************  «*«>*>  *«**4l*****:4t4[«**:»4l:4l* 


A.  15  therme.for 

c — get  thermodynamic  functions  (h,c_p,s-sO)  for  electrons. 

C  t  Itcitc  %  1)1  :|c  ***  %**>)(*  iK  It!  41 ***  ><1  IK**  4=  **** ’ll  ***  **  it:  ***#%*  * 

SUBROUTINE  thermeC 
ttemp .pressl , enthe , Cpe , entre) 

REAL*8  temp , press 1 , enthe , Cpe , entre 
IMPLICIT  REAL*8  (a-h.l-z) 

DATA  gasR, atiii/8 . 31439 , 0 . 101325e6/ 

DATA  tmpO  /2.9815E2/ 
enthe=2 . 5*gasR*temp 
Cpe=2.5*gasR 

entre=gasR*(2.5*DL0G(temp/tmp0)-DL0G(pressl/atin)) 

END 

C*****+it'>l'****'K*4i**i)i***4'***4i4i*>|i*>)i*4ii|ii|ii)ii|:i)i*i|i**>|t*i|t4i4i4ii)t4t4nt'4ti|ti)i*4i*4t4ti|ii|i4i*i):4t*4t** 


A. 16  rhsl.for 


c--gets  the  right-hand-sides  (rhs)  of  the  Saha  equation  for  use  in  the 
c — nuclear  fraction  (y)  calculations.  See  paper  theory  for  derivations, 
c — Say  there  aire  the  neutral  atomic  species,  4  ionized  atomic  species, 
c — and  the  two  diatomic  species:  7  species  per  element.  We  only  need  6 
c — equations,  and  therefore  6  rhs’s,  to  solve  for  these  7  species, 
c — since  the  seventh  equation  is  conservation  of  mass.  Therefore,  for 
c — the  7  species  we  have,  rhs (5)  will  be  left  empty  (0),  and  rhs (6)  and 
c — rhs(7)  correspond  to  the  diatomic  species  equations.  If  we  only  go 
c — up  to  A+3,  then  rhs (4)  and  rhs (5)  will  be  empty. 

C>|c:|c4c*i|c%4c  4c  4c4c4E**iti*iti*4i****iti4‘*itE**4iit‘**:4‘*4‘4[4t*******%*****!ti*****:tc>tc*itE**:|c*  **:**** 

SUBROUTINE  rhsK 

fttemp,pressl,ist,ion,Z,uD,uI,MW,rhs) 

INTEGER  ist 

REAL*8  temp,pressl,ion(5) ,Z(7,3) ,uD,uI.MW,rhs(7) 

INTEGER  i 

IMPLICIT  REALMS  (a-h,l-z) 

DATA  coefD,coefI.K_eV  /2. 59466E3, 3. 33381E-2. 11604.85/ 
tempinv=l . 4387/temp 
tempre=temp*temp*DSQRT(temp)/pressl 
DO  100  i=l,7 
100  rhs(i)=0. 

DO  200  i=l,ist-l 
boltz=0. 

IF  (ion(i)*tempinv.LT.70.)  boltz=DEXP(-ion(i)*tempinv) 
rhs(i)=coef I*Z(i+l,l)/ZCi,l)*boltz*tempre 
200  CONTINUE 

mu=MW*MW/(MW+MW) 

boltz=0, 

IF  (uD*K_eV/temp.LT.70.)  boltz=DEXP(-uD*K_eV/temp) 
rhs ( 6) =coef D*mu*DSQRT (mu) *boltz*tempr e* 
ftZ(l,l)t‘Z(l,l)/ 

4Z(6,1)/Z(6,2)/Z(6,3) 

bolt2=0. 

IF  (uI*K_eV/temp.LT.70. )  boltz=DEXP(-uI*K_eV/temp) 
rhs(7)=coef I*boltz*tempre* 
ftZ(7,l)*Z(7,2)*Z(7,3)/ 
fcZ(6,l)/Z(6,2)/Z(6,3) 

END  ! subroutine  rhsl 

C  4c  4c  itc  4: 4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  i|c  4c  *  ’ll  *  *  *  4c  4  4c  4c  4c  4c  ]|c  4c  itc  4citc  4c  4c  4c  4c  4  4  4c  4  4c  4  4c  4c  4  4c  4c  4c  4c  4c  4c  4citc  4c  etc*  4c  4c  4c  4c  4c  4c  4c 
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A. 17  rhs2.for 


c — gets  right-hand-sides  of  Saha  equations  for  diatomic,  heteronuclear 
c — species.  See  paper  theory  for  derivation. 

Q^ii/iin:^7lfmiti^m/iiHf^iif}tf*^fitf***************************************************** 

SUBROUTINE  rhs2( 

£temp , pressl . ZA . ZB , ZAB . RAB , lAB ,MWA , MWB , rhsAB) 

REAL+8  temp, pressl, ZA (7, 3) ,ZB(7,3) ,ZAB(2,3) ,RAB, lAB, 
4MWA,MWB,rhsAB(2) 

IMPLICIT  REAL*8  (a-h,l-2) 

DATA  coefD,coefI,K_eV  /2.59466E3.3.33381E-2,11604.85/ 

tempre=temp*temp*DSqRT(temp) /pressl 

boltz=0. 

IF  (RAB*K_eV/temp.LT.70. )  bolt2=DEXP(-RAB*K_eV/temp) 
muAB=MWA*MWB/ (MWA+MWB) 

rhsAB ( l)=coef D*muAB*DSQRT (muAB) *boltz*tempre* 
m(l,l)*ZB(l,l)/ 
ftZAB (1,1) /ZAB (1,2) /ZAB (1,3) 
boltz=0. 

IF  (IAB*K_eV/temp.LT.70.)  bolt2=DEXP(-IAB*K_eV/temp) 
rhsAB (2)=coefI*bolt2*tempre* 
ftZAB(2,l)*ZAB(2,2)*ZAB(2,3)/ 
ftZAB(l,l)/ZAB(l,2)/ZAB(l,3) 

END  ! subroutine  rhs2 

C ** itc 1)1  4c :tc %:)( itc If  Itc # itc itc  4cit( iK :*<**********  * :4c itc *  :4c *  3|c  *  :4c  4c  :tc « >tc Itf  :4c  Itc « 
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A. 18  gesl.for 

c — guesses  (hence  ‘‘ges”)  the  initial  species  distribution  for  single 
c — (hence  "1”)  element,  based  on  the  assumption  that  only  two  species 
c — dominate  at  at  T,P.  This  assumption  linearizes  the  otherwise  very 
c — nonlinear:  equations,  and  gives  the  right  answer  to  within  5*/,  for  a 
c — single  element,  or  for  an  inert  (no  chemical  reactions)  mixture  of 
c — multiple  elements.  Lickity  split. 

c — The  conditions  for  checking  which  two  species  are  dominant  are 
c — these:  say  we’re  examining  N_+l,  and  N_+3  (abbrv.  Nl,  N2,  N3) 

c — The  switchover  point  from  considering  N1,N2  to  considering  N2,N3 
c — comes  when  there  are  equal  concentrations  of  Nl  and  N3.  Using  this, 
c — the  fact  than  Nl+N2+N3=N_heavy,  and  the  fact  that 
c — Nl+2N2+3N3=M_electrons  all  give  simple  relations  for  the  crossover 
c — points;  what  makes  the  crossover  points  unique  for  an  element  are 
c — the  values  of  the  right-hand-sides  of  the  Saha  equations,  which 
c — depend  on  the  molecular  weight,  partition  functions,  etc. 

SUBROUTINE  gesl( 
trhs,yO,y) 

REAL*8  rhs(7),y0,y(7) 

IMPLICIT  REALMS  (a-h,l-z) 

IF  (rhs(6)*rhs(l).LT.le-3)  THEN 
y(l)=yO*DSQRT(rhs(6)/(4.+rhs(6))) 
y(6)=y0-y(l) 

ELSEIF  (rhs(l)*rhs(2).LT.0.25)  THEN 
y(2)=yO/DSqRT(l+rhs(l)) 
y(l)=y0-y(2) 

ELSEIF  (rhs(2)*rhs(3).LT.0.44)  THEN 
y(3)=y0*(-0.5+DSQRT(0.25+2.*rhs(2)/(l.+rhs(2)))) 
y(2)=y0-y(3) 

ELSEIF  (rhs(2)*rhs(3).GE.0.44)  THEN 
y(4)=yO*(-l.+DSqRT(l.+3.*rhs(3)/(l.+rhs(3)))) 
y(3)=y0-y(4) 

ENDIF 

CALL  getye(y,ye) 

IF  (ye.Eq.O.)  THEN 

y(2)=y(l)*rhs(l)+DSqRT(y(l)*y(l)*rhs(l)*rhs(l)+ 
ft  y(l)*rhs(l)*(y(l)+y(6))) 

y(7)=y(6)*rhs(7)+DSqRT(y(6)+y(6)*rhs(7)*rhs(7)+ 
ft  y(6)*rhs(7)+(y(6)+y(l))) 

ENDIF 

END  ! subroutine  gesl 


G6 


Q***^i^i**:it^**************************************************^^*********** 
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A. 19  itery.for 

c  find  which  nuclear  fraction  is  largest,  and  set  the  equations  to 
c — “pivot”  around  that  species. 

C*****************l(l********+]|l*%*)(t**************i(t********:(,:4,3(j^3^*^^3),,^j,j^jj,j^^ 

SUBROUTINE  iteryC 
ky.ity) 

INTEGER  ity 
REAL+8  y(7) 

IMPLICIT  REAL*8  (a-h.l-z) 
ity=6 

IF  (y(l).GT.y(6))  THEN 
ity=l 

ELSEIF  (y(2) .GT.y(l))  THEN 
ity=2 

ELSEIF  (y(3).GT.y(2))  THEN 
ity=3 

ELSEIF  (y(4).GT.y(3))  THEN 
ity=4 
END  IF 

END  (subroutine  itery 
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^  A. 20  getyh.for 

c — get  the  sum  of  the  nuclear  fractions=the  “heavy’’  fraction,  yh. 
C********************************************************************** 
SUBROUTINE  getyhC 
fty.yh) 

REAL*8  y(7).yh 
INTEGER  i 

IMPLICIT  REAL*8  (a-h.l-z) 
yh=0. 

DO  100  i=l,7 
100  yh=yh+y(i) 

END  {subroutine  getyh 

C*********************************^l4i**^i*^f*^i**:^ilf:t^t********************** 
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A. 21  getye.for 

c  get  the  electron  fraction,  ye,  based  on  the  charge  weighted  average 
c — fraction 

SUBROUTINE  getyeC 
fty.ye) 

REAL*8  y(7),ye 
INTEGER  i 

IMPLICIT  REAL*8  (a-h,l>z) 
ye=0. 

DO  100  i=l,5 

100  ye=ye+(i-l)*y(i) 
ye=ye+y(7) 

END  ! subroutine  getye 

C***************************^*****^*:^*-i.4.^f^f^,^i^)Ht**m********************* 
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^  A.  22  getwyh.for 

C"get  the  mass-weijhted  sum  of  the  heavy  nuclear  fractions.  Used  for 
c — conservation  of  mass  equations. 

SUBROUTINE  getwyhC 
&y,yAB,wyh) 

REAL+8  y(7).yAB(2).wyh 
INTEGER  i 

IMPLICIT  REAL+8  (a-h.l-z) 
wyh=0 . 

DO  100  i=1.5 
100  wyh=wyh+y(i) 

wyh=wyh+2 . *y (6) +2 . ♦y (7 ) +yAB ( 1 ) +yAB (2) 

END  ! subroutine  getuyh 

C********************************************************************** 
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A. 23  balance.for 


c — the  inner  loop.  Finishes  when  either  (1)  the  mass  weighted 
c — nuclear  fraction  sum  (vyh)  is  within  epsilon  of  the  total  nuclear 
c — fraction,  yO,  or  (2)  the  sum  of  the  heavy  fractions  (yhh,  a 
c — different  name,  so  we  don’t  screw  up  the  outer  loop)  is  very  very 
c — small.  This  protects  us  against  the  case  where  all  of  an  element 
c — combines  with  the  other  to  form  the  compound  AB;  otherwise,  the  routine 
c — will  keep  trying  to  reduce  the  partial  residual  by  making  the 
c — concentrations  smaller  and  smaller,  until  finally  we  get  an 
c — underflow  error. 

SUBROUTINE  balance ( 

*ye . yh . beta , epsilon , chi , ymin , yAB , seein , yO , y , ity , rhs) 

INCLUDE  format. for 
INTEGER  ity, seein 

REALMS  ye, yh, beta, epsilon, chi, ymin, yAB(2) ,y0,y(7) ,rhs(7) 

INTEGER  iin 

IMPLICIT  REAL*8  (a-h.l-z) 

loy=0. 

hiy=yO 

iin=0 

wyh=-l. 


yhh=l . 

c - the  inner  loop. 

DO  WHILE  ((DABS(wyh-yO).GT. epsilon). AND. (yhh. GT. ymin)) 

c - get  the  new  (“nu”)  fractions  (“y”) 

CALL  nuy(ity,yO,y,ye,yh,beta,rhs) 

c - then  get  the  new  weighted  heavy  fraction 

CALL  getwyh(y,yAB,wyh) 
c - under-relax/bisection  the  pivot  y 


IF  (wyh.GT.yO)  THEN 
hiy=y(ity) 

y(ity)=(l.-chi)*y(ity)+chi*loy 
ELSEIF  (wyh.LT.yO)  THEN 
loy=y(ity) 

y(ity)=(l .-chi)*y(ity)+chi*hiy 
ENDIF 


c - get  the  new  weighted  heavy  fraction  with  the  updated  pivot  y 

CALL  getwyh(y,yAB,wyh) 

c - and  get  the  test-case  heavy  fraction 

CALL  getyh(y,yhh) 
iin=iin+l 


c 


how’s  the  progress? 

IF  (mod(iin,seein).EQ.O)  THEN 
WRITE(*.192)’+’,iin.wyh-yO 
ENDIF 
ENOOO 
END 

C******************************* **************************** *********** 
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c — given  the  pivot  y  (ity  is  the  array  index  for  the  y  array),  the 
c — total  nuclecu:  fraction  yO,  the  temporary  ye, yh, beta,  and  the  rhs 
c — array,  we  update  the  other  nuclear  fractions  for  an  element. 
C*********************************************************!^* *********** 
SUBROUTINE  nuy( 
tity,yO,y,ye,yh,beta,rhs) 

INTEGER  ity 

REAL*8  y0,y(7) ,ye,yh,beta,rhs(7) 

IMPLICIT  REAL*8  (a-h.l-z) 

c - check  if  electron  density  is  zero,  and  should  still  be  zero.  Need 

c - to  check,  because  otherwise,  will  try  to  divide  by  beta,  which  is 

c - zero  if  ye=0. 

IF  (ye.EQ.O.)  THEN 
y(2)=0. 
y(3)=0. 
y(4)=o. 
y(6)=o. 
y(7)=o. 

c - calculate  the  other  non-charged  species  from  the  pivot  species 

IF  (ity.EQ.6)  THEN 
y(l)=DSQRT(rhs(6)*y(6)*(yO-y(6))) 

ELSEIF  (ity.EQ.l)  THEN 
y(6)=2.*yCl)*y(l)/rhs(6)/(yO+y(l)) 

ENDIF 

c - find  the  other  nuclear  fraction  y’s  from  the  pivot  y.  The 

c - equations  are  only  slightly  non-linear  when  we  (temporarily)  hold 

c - the  electron  fraction,  ye,  and  beta  (as  well  as  yh)  constant. 

ELSEIF  (ye.GT.O.)  THEN 
IF  (ity.EQ.e)  THEN 
y ( 1 ) =DSQRT ( ( yh+y e) *rhs (6)*y(6)) 
y(2)=rhs(l)/beta*y(l) 
y(3)=rhs(2)/beta*y(2) 
y(4)=rhs(.3)/beta*y(3) 
y(5)=rhs(4)/beta*y(4) 
y(7)=rhs(7)/beta*y(6) 

ELSEIF  (ity.EQ.l)  THEN 
y(2)=y(l)*rhs(l)/beta 
y(3)=y(2)*rhs(2)/beta 
y(4)=y(3)>t‘rhs(3)/beta 
y(5)=y(4)*rhs(4)/beta 
y(6)=y(l)*y(l)/rhs(6)/(ye+yh) 


y(7)=y(6)*rhs(7)/beta 
ELSEIF  (ity.EQ.2)  THEN 
y(l)=y(2)*beta/rhs(l) 
y(3)*y(2)*rhs(2)/beta 
y(4)=y(3)*rhs(3)/beta 
y(5)=y(4)*rhs(4)/beta 
y(6)=y(l)*y(l)/rhs(6)/(ye+yh) 
y(7)=y(6)*rhs(7)/beta 
ELSEIF  (ity.EQ.3)  THEN 
y(2)=y(3)*beta/rhs(2) 
y(l)=y(2)*beta/rhs(l) 
y(4)=y(3)*rhs(3)/beta 
y(5)=y(4)*rhs(4)/beta 
y(6)=y(l)*y(l)/rhs(6)/(ye+yh) 
y(7)=y(6)*rhs(7)/beta 
ELSEIF  (ity.EQ.4)  THEN 
y(3)=y(4)*beta/rhs(3) 
y(2)=y(3)*beta/rhs(2) 
y(l)*y(2)*beta/rhs(l) 
y(5)*y(4)*rhs(4)/beta 
y(6)=y<l)*y(l)/rhs(6)/(ye+yh) 
y(7)*y(6)*rhs(7)/beta 


ENDIF 

ENDIF 

END  ! subroutine  nuy 

C********************************************ltli**^i*!ti****:tL**:lf:ff:4:^i>ti^ii^^^i}ti^i^c* 
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A. 25  resl.for 

c — get  the  residual  for  an  element,  subroutine  can  be  used  for  either 
c — element;  add  residuals  from  both  elements,  and  hetero-diatomics  to 
c — get  total  residual. 

Q>)i***4t**4i**4>iti4i:ti****i4i4i*****4i***4i*****>ti*****>t‘*«***  ************  *********** 

SUBROUTINE  resl( 

ftyO , y , rhs , yAB , ye , yh , beta , residual) 

REAL*8  y0,y(7) ,rhs(7) .yAB(2) , ye, yh, beta, residual 
INTEGER  i 

IMPLICIT  REAL*8  (a-h,l-z) 

RE1AL*8  res  (7) 

CALL  getwyhCy.yAB.wyh) 
res(l)=  wyh-yO 

res(2)=-rhs(l)*y(l)+beta*y(2) 

res(3)=-rhs(2)*y(2)+beta*y(3) 

res(4)=-rhs(3)*y(3)+beta*y(4) 

res(5)=-rhs(4)*y(4)+beta*y(5) 

res(6)=  y(l)*y(l)/(yh+ye)-rhs(6)*y(6) 

res(7)=-rhs(7)*y(6)+beta*y(7) 

residual=0 . 

DO  100  i=l,7 

residual=residual+DABS(res(i)) 

100  CONTINUE 

END  .‘subroutine  resl 

(^iti*****************************************  *****************  *********** 


A. 26  res2.for 

c — get  residual  for  heteronuclear  diatomic  species. 

SUBROUTINE  res2( 

&yA , yB , yAB , rhsAB , ye , yh , beta , residual) 

REAL*8  yA(7) ,yB(7) ,yAB(2) ,rhsAB(2) , ye, yh, beta, residual 
INTEGER  i 

IMPLICIT  REAL*8  (a-h,l-z) 

REAL*8  res (2) 

res(l)=  yA(l)*yB(l)/(yh+ye)-rhsAB(l)*yAB(l) 

res(2)=-rhsAB(2)*yAB(l)+beta*yAB(2) 

residual=0. 

DO  100  i=l,2 

res idual=res idual+DABS (res(i)) 

100  CONTINUE 

END  ! subroutine  res2 

Q^:ilf^f*^ilf:lf**:t:*:ltililf*ilt*******:t^i^*^f*******t*********************************** 
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A. 27  rocket. for 


c — get  performaince  from  an  ideal  rocket,  assuming  (1)  our 
c — temperature  and  pressure  are  that  of  the  rocket  chamber,  (2)  frozen 
c — flow  (no  change  in  equilibrium  species  composition  downstream),  (3) 
c — isentropic  (reversible)  expansion  of  gases,  (4)  hard  sphere  ideal  gas 
c — (so  forget  all  the  complicated  thermodynauaic  functions  we  calculated; 
c — just  use  c_p=2.5R).  Frozen  flow  assumption  relieves  us  from  having 
c — to  calculate  coronal  equilbrium  at  exit  conditions. 

SUBROUTINE  rocket ( 

4yOA,yOB,MWA,MWB,ye,yh,temp,pressl.press2,mdot,thrust , Isp,Pjet) 
REAL+8  yOA,yOB,MWA,MWB,ye,yh, 

At emp,pressl,press2,mdot, thrust, Isp,P jet 
REAL*8  MWbar, prat, vex 
DATA  gasR.grav/8. 31439, 9. 8062/ 

MWbar=(yOA*MWA+yOB*MWB)/(yh+ye) 

prat=press2/pressl 

vex=DSQRT (5 . e3*gasR/MWbar*temp* ( 1-prat *prat  4>DSQRT (prat ) ) ) 

thrust=mdot*vex 

Isp=vex/grav 

P j  et=0 . 5*mdot*vex*vex 

END  .'subroutine  rocket 

C****4<4t4>4<*4<*4>4i4<**4c*4i4<4<4<*4c*4i4c*4c4i4:4i*4c**4<4c4c4c4‘4c4c4c4<4c4c4:4c4c4c44<4c4c4c4‘4<4c4c4c44c4<4c4<4c4<4c4c 
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A. 28  output.for 

c — The  mother  of  all  subroutines.  Basically  has  two  branches:  (1)  if 
c — the  has  selected  howmany=l,  this  subroutine  will  print  anything  of 
c — interest  to  a  file  name  name. out,  where  "name''  is  the  name  of  the 
c — compound  formed  by  the  two  species  the  user  chose.  The  user  gets  his 
c — original  input  data,  the  electronic,  rotational,  and  vibrational 
c — partition  functions  of  all  the  species,  the  right  hand  sides  of  the 
c — saha  equation,  the  nuclear  fractions  (y's)  of  the  species,  the  mole 
c — fractions  of  the  species,  the  densities,  the  heats  of  formation  for 
c — the  species,  the  enthalpies,  heat  capacities,  and  entropies,  both  of 
c — the  species  emd  for  the  entire  gas,  and  the  rocket  performance.  If 
c — you  want  it  and  its  not  there,  this  program  doesn't  do  it.  (2)  If  the 
c — user  selects  multiple  temperatures  (howmany.!=  1),  a  more 
c — conventient,  tabulated  form  of  data  is  outputted  to  several  files, 
c — one  file  per  data  type  per  element.  However,  due  to  the  MS-DOS 
c — limitations  on  the  number  of  possible  files  open  at  one  time,  the 
c — user  must  choose  between  density  and  mole  fraction.  I  chose  these 
c — since  it’s  not  likely  the  user  will  want  both.  Depending  on  the 
c — user’s  choice  of  howmany,  subSUBRQUTINEs  open  the  appropriate  files 
c — and  make  the  appropriate  headers. 

QilHi/:if^l*!ti^l7tt*if*1fJtt*t******************************if************************ 

SUBROUTINE  output ( 

tyOA , yOB , tempi , temp2 , howmany , pressl , press2 , mdot , temp , itemp , 

&chA , DA2 , IA2 , ionA , ZA , rhsA , yA , enthA , CpA , entr A , 

4chB , DB2 , IB2 , ionB , ZB , rhsB , yB , enthB , CpB , entrB , 

ftname , RAB , lAB , ZAB , rhsAB , yAB , enthAB , CpAB , ent rAB , enthe , Cpe , entre , 

&Pjet, thrust, Isp,densmol) 

CHARACTER  chAx‘2,chB*2,name*4 
INTEGER  itemp, densmol, howmany 
REAL*8 

ftyOA , yOB .tempi , temp2 , pressl , press2 ,mdot , temp , 

fcDA2,IA2,ionA(5) .ZA(7,3) ,rhsA(7) ,yA(7) . enthA(7) ,CpA(7) , entrA(7) , 
ftDB2,IB2,ionB(5) .ZB(7,3) ,rhsB(7) .yB(7) .enthB(7) .CpB(7) .entrB(7) . 
tRAB,IAB,ZAB(2,3),rhsAB(2),yAB(2),enthAB(2) ,CpAB(2) .entrAB(2) , 
ftenthe , Cpe , entre , 
tP jet, thrust, Isp 
CHARACTER  filename*12 
INTEGER  i 

IMPLICIT  REALMS  (a-h,l-2) 

REALMS  dirA(7) ,dirB(7) ,dirAB(2) ,mf A(7) ,mfB(7) ,mf AB(2) , 
tdenA(7) ,denB(7) ,denAB(2) 

DATA  bolt2k/1.38066E-23/ 
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DATA  JmoleV, Jmolcni/96484 . 52. 11 . 96263/ 

INCLUDE  format. for 

re-evaluate  the  heavy  and  electron  fractions,  just  to  be  safe. 
CALL  getye(yA,yeA) 

CALL  getye(yB,yeB) 

CALL  getyh(yA,yhA) 

CALL  getyhCyB.yhB) 
ye=yeA+yeB+yAB (2) 
y h=y hA+y hB+y AB ( 1 ) +yAB ( 2 ) 
beta=ye/(yh+ye) 

pkt=pressl/(boltzk*temp*(ye+yh)) 
enthal=0 . 
heatcp=0. 
entrop=0 . 

calculate  the  heats  of  formation  for  the  species: 

dissociation  (“d"),  ionization  (“i”),  and  reactions  (“r”) 

also  calculate  the  gas  thermodynamic  properties,  based  on  the 

species  properties  and  the  mole  fractions. 

dirA(l)*0 . 5*DA2* JmoleV 

dirB(l)=0.5*DB2* JmoleV 

DO  100  i=l,5 

enthal=enthal+yA ( i ) *enthA ( i ) +yB ( i ) ♦ent hB ( i ) 
heat cp=heat  cp+yA ( i ) *CpA ( i ) +yB ( i ) *CpB ( i ) 
entrop*entrop+yA(i)*entrA(i)+yB(i)*entrB(i) 
dirA(i+l)=dirA(i)+ionA(i)*Jmolcm 
dirB(i+l)=dirB(i)+ionB(i)*Jmolcm 
mfA(i)*yA(i)/(yh+ye) 
mfB(i)=yB(i)/(yh+ye) 
denA (i)=yA(i) ♦pkt 
denB ( i ) =yB ( i ) *pkt 
CONTINUE 
DO  200  i=l,2 
enthal=enthal+ 

k  yAB(i)*enthAB(i)+yA(i+5)*enthA(i+5)+yB(i+5)*enthB(i+5) 
heatcp=heatcp+ 

k  yAB(i)*CpAB(i)+yA(i+5)*CpA(i+5)+yB(i+5)*CpB(i+5) 
entrop=entrop+ 

k  yAB(i)*entrAB(i)+yA(i+5)*entrA(i+5)+yB(i+5)*entrB(i+5) 
mfA(i+5)=yA(i+5)/(yh+ye) 
mfB(i+5)=yB(i+5)/(yh+ye) 
denA ( i+5) =yA ( i+5) *pkt 
denB ( i+5) =yB (i+5) *pkt 


mf AB(i)=yAB(i)/(yh+ye) 
denAB ( i ) =yAB ( i ) *pkt 
CONTINUE 
dirA(6)=0. 
dirB(6)=0. 
dirA(7)=IA2*JmoleV 
dirB(7)=IB2*JmoleV 
dirAB(l)=(0.5*(DA2+DB2)-RAB)*JmoleV 
dirAB(2)=(0.5*(DA2+DB2)-RAB+IAB)*JmoleV 
eiithal=enthal+ye*enthe 
heatcp=heatcp+ye*Cpe 
entrop=entrop+ye*entre 
enthal=eiithal/(ye+yh) 
heatcp=heatcp/ (ye+yh) 
entrop=entrop/ (ye+yh) 

-send  a  teaser  to  the  screen 
WRITE(*.191)'  Pjet  (W)  ’.Pjet 
WRITEC*, 192) ’  Thrust (N) , Isp(s) * .thrust, Isp 
I F  ( howmany . EQ . 1 )  THEN 
IF  (itemp.EQ.l)  THEN 

OPEN  (UNIT=90.  FILE=filename (name, '.out')) 

ENDIF 

WRITE(90 , ' (A) ’ ) ' Using  ' //chA// ' , ' //chB// ' , ’ //name 
WRITE(90,*) 'Stoichiometric  Proportions  of  '//chA// ' , ’//chB 
WRITE(90.’("  [”.1PE11.3,".",1PE11.3,"]  ")’)yOA.yOB 
WRITE(90,*) 'chamber,  exit  pressures  (atm)' 

WRITE(90,'("  [",1PE11.3.  ",",1PE11.3,  "]  " )  ’ )pressl .press2 
WRITE(90,*) 'Temperature  range:  tempi (K) ,temp2(K) .howmany ' 
WRITE(90,'("  [".1PE11.3.",",1PE11.3,",",I5,"]  ")') 
ft  tempi, temp2, howmany 

WRITE(90,*) 'mass  flow  rate  (kg/s)’ 

WRITE(90,’(”  [",1PE11.3,’']  ”)’)mdot 

- pretty  display  of  partition  function 

CALL  Zlshow(chA,ZA) 

CALL  Zlshow(chB,ZB) 

CALL  Z2show(name,ZAB) 

- show  the  information  of  a  species,  described  by  the  first 

- parameter  string 

CALL  specshow( 'Right-hand-sides  of  Saha  Equation’, 
ft  chA,chB,n2une,rhsA,rhsB,rhsAB) 

CALL  specshow( 'Nuclear  fractions.y' , 
ft  chA, chB, name, yA,yB,yAB) 


WRITEOO ,  194)  * yh , ye , yh+ye , beta ’ , yh .  ye ,  yh+ye , beta 
CALL  specshowC’Mole  fractions, y/(yh+ye) ’ , 
k  chA.chB.name.mf A.mfB.mf AB) 

WRITEOO, 191)  ’e-’  ,ye/(yh+ye) 

CALL  specshowC ’Species  densities  (m'-S) ’ , 
k  chA,chB,name,denA,denB,denAB) 

WRITEOO,  191)  ’e-’  .pkt*ye 

CALL  specshowC’Heat  of  formation  (J/mol) ’ , 

4  chA,chB,name,dirA,dirB,dirAB) 

CALL  specshowC ’Enthalpies  (J/mol)’, 
k  chA,chB,name,enthA,enthB .enthAB) 
WRITE(90,191)’e-’,enthe 
WRITE(90,191)name//’  gas’,enthal 
CALL  specshowC ’Heat  capacities  (C_p)  (J/K/mol)’, 
k  chA,chB,name,CpA,CpB,CpAB) 

WRITE(90,191)’e-’,Cpe 
WRITE(90,191)name//’  gas’,heatcp 
CALL  specshowC ’Entropy  diff  (from  STP)  (J/K/mol)’. 
k  chA,chB,name,CpA,CpB,CpAB) 

WRITE(90,191)’e-’.entre 
WRITE(90,191)name//’  gas’,entrop 
WRITE(90.*) 

WRITE(90,192) ’Thrust (N) ,Isp(s) ’ , thrust, Isp 
WRITE(90.191)’Pjet  (W)’,Pjet 
ELSE 


IF  (itemp.EQ.l)  THEN 

c - open  all  the  appropriate  files  for  output 

CALL  opunhdr ( chA , chB , name , densmol ) 

ENDIF 

c - write  all  the  data  to  all  the  files 

WRITE(61 , 139) temp , enthA , enthe 


WRITE (62,139 ) t emp , enthB , enthe 
WRITE(63, 134)temp, enthAB, enthe 
WRITE (71,139)t  emp , CpA , Cpe 
WRITE(72,139)temp,CpB,Cpe 
WRITE(73 , 134) temp , CpAB , Cpe 
WRITE(81 , 139) temp , entrA , entre 
WRITE(82 , 139) temp , entrB , entre 
WRITE(83 , 134) temp , entrAB , entre 

c - the  MS-DOS-imposed  choice 

IF  (densmol. EQ.l)  THEN 
WRITE(91, 139) temp, denA,pkt*ye 
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WRITE(92, 139) temp, denB,pkt*ye 
WRITE(93 , 134) temp , denAB , pkt*ye 
ELSE 

WRITE(91 , 139) temp ,mf A , ye/ (yh+ye) 

WRITE(92, 139)temp,mfB.ye/(yh+ye) 

WRITE ( 93 , 134) t  emp , mf  AB , ye/ (yh+ye) 

ENDIF 

WRITE(94 , 134) temp , P j  et , thrust , Isp 
ENDIF 

END  ! subroutine  output 

C4c]tc**4E4i**#*4i:<:***4<***!ti4:*4!***4i***4:*^*«4i*«*«*******  ********’)<* ’t‘***’*‘****  +  +  +  * 
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A. 29  opunhdr.for 

c — opens  (“opun”)  all  the  output  files  for  the  howmany!=l  case,  and 
c — writes  the  headers  (“hdr"),  such  as  a  line  describing  what  the  file 
c — contains,  with  units,  and  species  column  headings, 
c — file  unit  key 
c — 60’s=€nthalpy 
c — 70’s=heat  capcity 
c — 80 ’ s=entropy 
c — 90 ’ s=density/molef raction 
c — +l=element  A 
c — +2=element  B 
c — +3=hetero-dia  molecule  AB 
c — output  file  naming  convention  is  as  follows: 

C"(name)_(description)  .out 

c — where  (name)  is  the  element  or  hetero-dia  molecule  name 
C”and  (description)  is  ■(eth,cp,etr,den,mol}  for  enthalpy,  etc. 

Qif^f**:^^**:^:^^^::^^*********************************************************** 

SUBROUTINE  opunhdr( 
tchA , chB , name , densmol) 

CHARACTER  chA*2,chB*2,name*4, 
ftf ilename*12,pref ix*4(3) ,exten*8(4) ,header*40(4) 

INTEGER  densmol.i, j ,ijunit 

pref ix(l)=chA 

pref ix(2)=chB 

pref ix(3)=name 

exten(l)='_eth.out’ 

exten (2) = ’ _Cp . out ’ 

exten(3)= ’ _etr . out ’ 

header (1)=’ Enthalpy  (J/mol)’ 

header(2)=’Heat  capacity,  C_p  (J/K/mol) ’ 

header(3)=’ Entropy  difference  from  STP  (J/K/mol)' 

IF  (densmol. EQ.l)  THEN 
exten(4)='_den.out ' 
headei (4)® 'Partial  density  (m“-3)' 

EL-E 

exten (4) = ' _mol . out ' 
header (4)=' Mole  fraction' 

ENDIF 

DO  10  i=l,4 
DO  10  j=l,3 

ijunit=60+10*(i-l)+j 

OPEN  (UNIT=ijunit,FILE=f ilename(pref ix(j) ,exten(i))) 
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10  CONTINUE 

OPEN  (UNIT=94,  FILE=f ilename (name, ’_rkt. out’)) 

DO  20  i=1.4 
DO  20  j=l,3 

ijunit=60+10*(i-l)+j 

20  WRITE(ijunit,*)pref ix(j) ,header(i) 

WRITE ( 94, ♦) neune,  ’  rocket  performance  parcimeters’ 

991  F0RMAT(2x,’temp(K)’,4x,A,9x,A,’+’,8x,A, ’+2’ ,7x,A, ’+3’ ,7x,A, 

ft'+4’,7x,A,’2’,8x,A,’2+’,7x,’e-’) 

993  F0RMAT(2x, ’tempCK) ’ ,4x, A,7x,A, ’+’ ,6x, ’e-’) 

DO  30  i=l,4 

ijunit=60+10*(i-l) 

WRITEC i j  unit+1 , 991 ) chA , chA , chA , chA , chA , chA , chA 
WRITE( i j  unit+2 , 991 ) chB , chB , chB , chB , chB , chB , chB 
WRITEC ijunit+3 , 993) name , name 
30  CONTINUE 

994  FORMAT (2x, ’temp(K) ’ ,4x, ’Pjet (W) ’ ,4x, ’ thrust (N) ’,2x,’Isp(s)’) 
WRITE(94,994) 

END  ! opunhdr 

C**************************************i^*ifit**************************** 
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A. 30  filename.for 


c — takes  two  strings  (of  arbitrary  length)  and  concatenates  them,  taking 
c — the  spaces  out  of  the  first  string.  This  is  particularly  useful  when 
c — making  a  filename  for  the  OPEN  command  that  MS-DOS  will  recognize, 
c — If  there  are  spaces  in  between  the  prefix  and  the  extension,  MS-DOS 
c — meOtes  a  file  with  no  extension... 

Q  #  4: Incite  :tc  :4l  411)1  %  Itc  1(1  Itc  4:  4l  Itt  4e  4ll|l  :tc  *  It!  Ik  Itr  iti  it  Ik  <|c  :|i  :4c  4c  ***  * 

CHARACTER*12  FUNCTION  f ilename(neime,exten) 

CHARACTER  name*(*) ,exten*(*) ,temp*12 

INTEGER  ilast 

temp=name 

ilast=INDEX(name,’  ') 

IF  (ilast. eq.O)  THEN 

temp(LEN(name)+l; )=exten 
ELSE 

temp(ilast : )=exten 
ENDIF 

f ilename=temp 
END 

Q^^^|^^^^|^Hi^Li)t*************************************************************** 
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A. 31  specshow.for 

c — displays  the  contents  of  a  species  variable  to  the  howinany=l  file, 
c--but  in  an  user-interpretable  fashion. 

SUBROUTINE  specshowC 
fcstr , chA , chB , name , yA , y B , yAB ) 

CHARACTER  str* (♦) , chA*2 , chB*2 , name*4 
REAL*8  yA(7),yB(7).yAB(2) 

100  FORMAT  (4x,5(’atomic, ’ ,4x) ,2(’diatomic, ’ ,2x)) 

101  FORMAT  (4x, ’neutral’ ,4x,4('+’ ,11, 9x) , ’neutral’ ,4x, ’+!’) 

102  FORMAT  (A.1P7(E11.3)) 

103  FORMAT  (A,53x,lP2(E11.3)) 

WRITEOO,*) 

WRITE(90,*)str 

WRITE(90,100) 

WRITE(90, 101)1,2,3.4 

WRITE(90,102)chA.yA 

WRITE(90,102)chB,yB 

WRITE(90,103)name,yAB 

END 

Qif:^*1f***^t1f^i*nn^***^i^**^t************************************************* 
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A. 32  Zlshow.for 


c  shows  the  peurtition  function  for  an  element  in  an  intelligible  way. 

SUBROUTINE  ZlshowC 
kch.Zl) 

CHARACTER  ch*2,partstr*3(3) 

REAL*8  Zl(7,3) 

INTEGER  i 

DATA  partstr/’el  ’ , ’vib’ , ’rot V 

100  FORMAT  (4z,5(’atomic, ’,4x),2(’diatomic, ’,2x)) 

101  FORMAT  (4x, ’neutral’,4x,4(’+’, II. 9x). ’neutral’, 4x. '+!’) 
WRITEOO.*) 

WRITEOO,*) ’Partition  Function:  ’,ch 
WRITEOO,  100) 

WRITEOO. 101)1, 2. 3, 4 
DO  10  i=l,3 

WRITEOO.  ’(A,1PE10.3,1P6(E11.3))’) 
ft  partstr(i).Zl(l,i),Zl(2,i),Zl(3.i), 
ft  Zl(4.i).Zl(5.i).Zl(6.i).Zl(7.i) 

10  CONTINUE 
END 

C*****l|<*****i(t*i(ii|!***:)t:(t***i^i(t**)(i^:4t************  ***!)[#  1(1*****  **j(t****i(c*****i(t*:(i* 
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A. 33  Z2show.for 


c — shows  the  hetero-diatomic  species  partition  function;  for  use  with 
c — the  howmany=l  output  file. 

C******************************* **♦******♦♦**♦♦♦*********♦****♦** ****** 
SUBROUTINE  Z2show( 

&name.Z2) 

CHARACTER  naune*4,partstr*3(3) 

REAU8  Z2(7.3) 

INTEGER  i 

DATA  partstr/’el  ’ , *vib’ , ’rot */ 

100  FORMAT  (4x,5(’atomic, ’ ,4x) ,2(’diatomic, ’ ,2x)) 

101  FORMAT  (4x, ’neutral’ ,4x,4(’+’ ,11, 9x) , ’neutral’ ,4x, ’+!’) 
WRITEOO,*) 

WRITE(90,*) ’Partition  Function:  ’.name 
WRITE(90,100) 

WRITE(90, 101)1.2.3.4 
DO  10  i=1.3 

WRITE(90, ’(A,54x.lP2(E11.3))’)partstr(i).Z2(l,i),Z2(2,i) 

10  CONTINUE 
END 

C********************************************************************** 
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A. 34  shutdown.for 


C’^’-shuts  down  the  program,  closing  all  the  open  files,  telling  the  user 
c-’-which  files  to  look  for,  and  waiting  for  an  <Enter>  before  returning 
c— control  to  DOS.  This  helps  when  running  the  prograun  from  a  DOS 
c— shell,  like  DCOM,  which  may  automatically  go  back  to  the  interface 
c—screen  and  prevent  the  user  from  seeing  which  files  he  needs  to  see. 

c  *  *  *  4c  ♦  lit  iti  ♦  ♦  ♦  *  If  *  If  ♦  *  *  Iti  *  4i  #  4^  ♦  *  *  m  >•<  *  *  *  *  itE  *  *  *  1*1  *  lie  *  *  *  *  %  4^  *  * 

SUBROUTINE  shutdown ( 
thowmany , chA , chB , name , densmol) 

CHARACTER  chA*2,chB*2,name*4, 
ftf ilename*12,pref ix*4(3) ,exten*8(4) ,pause*l 
INTEGER  densmol,i, j .howmany 
IF  (howmany.EQ.l)  THEM 

WRITEC*,*) ’Look  for  output  in  file  (’ ,f ilenameCname, ' .out') , ’ 
ELSE 

prefix(l)»chA 
pref ix(2)=chB 
prefix (3) =name 
exten ( 1 ) = ' _eth . out ’ 
exten(2)='_Cp.out' 
exten(3)='.etr.out’ 

IF  (densmol. EQ.l)  THEN 
exten(4)»’_den.out’ 

ELSE 

exten (4) = ’ .mol . out ’ 

ENDIF 

WRITEC*,*) 'Look  for  output  in  files:’ 

DO  10  i=l,4 
DO  10  j=l,3 

WRITEC*,*) ’(’, filename (pref ix(j), exten (i)) , ’) ' 

10  CONTINUE 

WRITE(*,*) ’ ( ’ , filename (name, ’_rkt .out ’) , ’) ’ 

ENDIF 

CL0SE(UNIT=61) 

CL0SE(UNIT=62) 

CL0SE(UNIT=63) 

CL0SE(UNIT»71) 

CL0SE(UNIT=72) 

CL0SE(UNIT=73) 

CL0SE(UNIT=81) 

CL0SE(UNIT=82) 

CL0SE(UNIT*83) 
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CL0SE(UNIT=90) 

CL0SE(UNIT=91) 

CL0SE(UNIT=92) 

CL0SE(UNIT=93) 

WRITEC*,*) ’Press  <Enter>  to  exit  program...’ 

READ(*, ’(A)’)pause 
END  [subroutine  shutdown 

C********************************************************************** 
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A. 35  format.for 


102 

FORMAT 

(214) 

111 

FORMAT 

(A. 14) 

112 

FORMAT 

(A. 214) 

113 

FORMAT 

(A. 314) 

131 

FORMAT 

(1P1(E11.3)) 

132 

FORMAT 

(1P2(E11.3)) 

133 

FORMAT 

(1P3(E11.3)) 

134 

FORMAT 

(1P4(E11.3)) 

135 

FORMAT 

(1P5(E11.3)) 

136 

FORMAT 

(1P6(E11.3)) 

137 

FORMAT 

(1P7(E11.3)) 

138 

FORMAT 

(1P8(E11.3)) 

139 

FORMAT 

(1P9(E11.3)) 

1310 

FORMAT 

(1P10(E11.3)) 

1311 

FORMAT 

(IPlKEll.S)) 

1312 

FORMAT 

(1P12(E11.3)) 

1313 

FORMAT 

(1P13(E11.3)) 

1314 

FORMAT 

(1P14(E11.3)) 

1315 

FORMAT 

(1P15(E11,3)) 

1316 

FORMAT 

(1P16(E11.3)) 

1317 

FORMAT 

(1P17(E11.3)) 

1318 

FORMAT 

(1P17(E11.3)) 

191 

FORMAT 

(A.lPEll.S) 

192 

FORMAT 

(A.1P2(E11.3)) 

193 

FORMAT 

(A.1P3(E11.3)) 

194 

FORMAT 

(A.1P4(E11.3)) 

195 

FORMAT 

(A.1P5(E11.3)) 

196 

FORMAT 

(A,1P6(E11,3)) 

197 

FORMAT 

(A.1P7(E11.3)) 

198 

FORMAT 

(A,1P8(E11.3)) 

199 

FORMAT 

(A,1P9(E11.3)) 

1910 

FORMAT 

(A,1P10(E11.3)) 

1911 

FORMAT 

(A,1P11(E11,3)) 

1912 

FORMAT 

(A,1P12(E11.3)) 

1913 

FORMAT 

(A,1P13(E11.3)) 

1914 

FORMAT 

(A,1P14(E11.3)) 

1915 

FORMAT 

(A.1P15(E11.3)) 

1916 

FORMAT 

(A,1P16(E11.3)) 
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9  B  Input  Files 


B.l  ABX.in 


ist  HWCaara) 

uD(eV) 

ul(eV) 

Tr.O(K) 

Tv,0(K) 

Tr.lCK) 

Tv.KK) 

c 

5 

12.011 

6.21 

12.1 

2.618 

2668. 

2.386 

1942. 

H 

2 

1.0079 

4.47 

15.4 

87.55 

6332. 

43.45 

3339. 

If 

4 

14.0067 

9.75 

15.5 

2.875 

3393. 

2.779 

3175. 

0 

8 

15.9994 

5.11 

12.0 

2.079 

2273. 

2.433 

2740. 

94 


B.2  ABXY.in 


uD(eV) 

ul(eV) 

Tr.OCK) 

Tv.O(K) 

Tr.KK) 

Tv.l(K) 

cc 

6.21 

12.1 

2.618 

2668. 

2.386 

1942 

CH 

3.46 

10.6 

20.80 

4112. 

20.39 

3941 

CN 

7.70 

14.1 

2.733 

2976. 

2.728 

2925 

CO 

11.0 

14.0 

2.778 

3121. 

2.844 

3185 

HH 

4.47 

15.4 

87.55 

6332. 

43.45 

3339 

HN 

3.47 

13.6 

24.02 

4722. 

22.08 

4204 

HO 

4.39 

12.9 

27.20 

5377. 

24.16 

4479 

NN 

9.75 

15.5 

2.875 

3393. 

2.779 

3175 

NO 

6.49 

9.26 

2.405 

2739. 

2.873 

3419 

00 

5.11 

12.0 

2.079 

2273. 

2.433 

2740 
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B.3  C_0.in 


90878.3 
1  0.0 
3  16.4 
5  43.5 
5  10193.70 
1  21648.4 
5  33735.2 
1  60333.80 
3  60353.00 
5  60393.52 
3  61982.20 
7  64088.56 
5  64093.19 
3  64092.01 
3  68858. 

3  69689.79 
5  69710.99 
7  69744.40 
3  70744.26 
1  71352.81 
3  71365.23 
5  71385.70 
5  72611.06 
1  73976.23 
9  75256.3 
5  77680.5 
1  78105.23 
3  78117.06 
5  78148.36 
5  78199.34 
7  78215.82 
9  78250.22 
3  78300.8 
5  78307. 

7  78316. 

3  78338. 

7  78531. 

3  78727.91 
5  79311.10 
3  79319.06 
1  79323.32 


3  80173.29 
5  80192.49 
7  80222.74 
3  80563.57 
3  81105.70 
1  81311.52 
3  81326.33 
5  81344.48 
5  81770.36 
1  82252.31 
5  83500. 

7  83761. 

3  83830. 

5  83837. 

7  83847. 

3  83882.5 
7  83949. 

3  84032. 

5  84102.6 
3  84112. 

3  84852.13 
5  84952. 

7  84986.2 
5  85400.38 
1  85625.84 
5  86187. 

5  86319. 

7  86326.7 
5  86371.3 
7  86396. 

3  86413.96 
7  86450. 

3  86491. 

5  86504. 

3  86517. 

5  87632. 

5  87706. 

7  87713. 

5  87752. 

7  87773. 

3  87795.3 
>  7  87807. 
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5  «7830  . 

3  87839. 

3  87631.3 
5  88541.8 
7  88547. 

7  88607. 

7  88624. 

3  88632.44 
5  88639. 

7  89081. 

5  89082. 

7  89146. 

7  89155. 

5  89158. 

5  89450. 

7  89514. 

7  89517. 

7  89779. 

7  89968.4 
-1  -1 


B.4  C_l.in 

196659.0 
2  0.0 
4  64.0 

1  46000.2 
4  43021.8 
6  43050.7 
6  74930.9 
4  74933.2 

2  96494.1 
2  110625.1 
4  110666.3 
2  116537.88 
2  131724.68 
4  131735.81 
4  142024.4 
4  145549.99 
6  145551.44 
6  150462.8 
4  150467.9 
2  157234.43 
2  162518.70 
4  162524.62 
2  166964.70 
4  166988.46 
6  167033.43 
4  168123.92 
6  168124.33 
2  168731.6 
4  168750.2 
14  168979.05 
2  173348.18 
2  175287.9 

4  175295.2 
2  178194.1 
4  178220.8 
6  178494.8 
14  178956.46 
2  181258. 

2  181694.50 
4  181709.20 
6  181734.21 
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8  181770.48 
2  182025.0 
4  182044.5 
6  184064.9 
14  184376.20 
4  184688.69 
2  186425.02 
4  186441.32 
6  186463.75 
4  188579.3 
6  188612.7 
2  194571.9 
4  195750.8 
6  195765.1 
8  195784.7 
10  195812.3 
2  196556.2 
4  196561.8 
6  196570.5 
8  196580.8 
-1  -1 


B.5  C-2.in 


386159.7 
1  0.0 
1  52315.0 
3  52338.0 
5  52394.8 
3  102351.4 
1  137374.0 
3  137403.4 
5  137450.5 
5  145875.1 
1  182520.2 
3  238160.7 
1  247169.5 
3  258931.4 
1  259653.8 
3  259659.3 
5  259672.1 
3  269957.6 
5  269959.7 
7  269962.9 
5  276482.7 
1  308162.9 
3  308196.2 
5  308264.8 
3  309404.5 

3  310005.2 
1  311720.7 

4  317743. 

5  317748. 

3  319719.4 
3  321358.8 
5  321375.1 
7  321398.6 
5  321949.1 
7  321955.8 
9  321964.7 
3  322403.1 
7  322701.1 
3  323024.0 
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5  323049.4 
7  323088.2 
5  324212.0 
3  327225.7 
1  329633.1 
3  329654.2 
5  329690.9 
5  332690.3 
5  333116.4 
5  333333.4 
7  333358.4 
9  333395.0 
3  337602.9 
5  337616.4 
7  337636.7 
3  339881. 

5  340049.5 
3  340075.8 
1  340090.3 
7  341368.5 
3  343255.7 
5  344181. 

I  345093.9 
7  345444. 

9  346525.1 

II  346253.0 
9  346577.5 
5  346656.0 
3  346713.1 
5  347099.5 
7  347101.3 
9  347103.7 
7  348859.5 
3  354796. 

3  357088. 

7  358046. 

9  358638.3 
11  358639.0 
9  358688.9 
5  358725.5 
9  358800. 
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7  359122.2 
3  363561. 

3  364896. 
15  365585. 
5  366027.0 
3  369926. 
15  370438. 
15  373748. 
5  376637. 

3  381104.8 
5  381919. 

7  381958. 

3  384313. 

5  384350. 

5  385637.5 
5  385816.2 
-1  -1 


B.6  C_3.in 


520177.8 
2  0.0 
2  64484.2 
4  64591.3 
2  302847.9 
2  320048.5 
4  320080.0 
4  324880.2 
6  324890.9 
2  401346.7 
2  408308.9 
4  408322.2 
4  410333.8 
6  410338.2 
8  410434.1 
2  445366.1 
2  448854. 

4  448861. 

6  449887.4 
14  449938.2 
18  449948.4 
2  468765. 

6  470763. 

10  471368. 
14  471403.0 
18  471407.4 
22  471407.9 
2  482659. 

6  483931. 

10  484309. 
14  484343.8 
18  484346.6 
22  484346.9 
6  492473. 

14  492743. 
36  492745. 
-1  -1 
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B.8  H.O.in 

109678.758 
2  0.000 
8  82259 
14  97492.3 
22  102823.9 
32  105291.6 
44  106632.15 
58  107440.42 
74  107965.03 
-1  -1 
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B.9  H-l.in 


100000000.000 
1  0.000 
-1  -1 
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B.10  N_0.m 

117345 
4  0.000 
6  19223 
4  19231 
6  28840 
2  83285.5 
4  83319.3 
6  83368 
2  86131.4 
4  86223.2 
6  88109.5 
4  88153.4 
2  88173 
2  93582.3 
2  94772.2 
4  94794.8 
6  94832.1 
8  94883.1 
2  95476.5 
4  95494.9 
6  95533.2 
4  96751.7 
4  96788.2 
6  96864.2 
2  97770.1 
4  97805.8 
6  99665 
4  99658 
2  103618.1 
4  103668.1 
6  103736.8 
2  104142.2 
4  104227.4 
4  104615.4 
2  104654.9 
4  104665 
6  104684 
8  104718 
10  104767 
6  104810.9 
8  104882.7 


108 


2  104864 
4  104890 
6  104957 
2  104987 
4  104998 
6  105011 
8  105020 
4  105120.8 
6  105144.3 
2  106478.6 
2  106760.5 
4  106780.1 
6  106816.1 
8  106870.7 
2  106982.7 
4  106998.3 
6  107039 
4  107447.2 
2  109813.5 
4  109857.8 
6  109927.9 
2  110029.2 
4  110108.5 
4  110196 
6  110214 
8  110248 
10  110304 
2  110221 
4  110275 
6  110288 
8  110339 
4  110221.7 
2  110244.6 
6  110311 
8  110373 
2  110325 
4  110351 
6  110403 
4  110448.3 
6  110470.5 
4  110521.9 
6  110545.8 


109 


2  112294.8 
4  112320.8 
2  112565.9 
4  112610.6 
6  112682.6 
2  112735 
4  112823 
4  112751 
6  112763 
8  112799 
10  112862 
4  112801 
2  112816 
6  112820 
8  112890.2 
6  112825 
6  112825 
8  112892 
2  112855 
4  112874 
6  112912 
4  112929.2 
6  112947.5 
2  114015 
4  114072 
6  114146 
2  114130 
4  114163 
28  114160 
12  114182 
8  114248 
4  114193 
2  114209 
6  114196 
8  114275 
4  114232.2 
6  114290.5 
6  114259 
6  114274 
2  114809 
4  114890 
6  114942 


110 


6  114950 
20  114988 
14  115004 
-1  -1 


111 


B.ll  N_l.in 


238846.7 
1  0.000 
3  49.1 
5  131.3 
5  15315.7 
1  32687.1 
5  47167.7 

7  92237.9 
5  92251.3 
3  92252.9 

8  109218.2 
1  109224.8 
5  144189.1 
1  148909.37 
3  148940.97 
5  149077.33 
3  149188.74 
3  155129.9 
3  164611.6 
3  166522.48 
5  166583.26 
7  166679.45 
3  166765.7 
3  168893.04 
1  170573.38 
3  170608.63 
5  170667 
5  174212.93 
1  178274.17 
5  186512.38 
7  186571.8 

9  186653.35 
5  187092.2 
3  187438.34 
5  187462.38 
7  187492.72 
5  188858.09 
3  188909.89 
1  188937.95 
7  189336 
3  190121.15 

i 


112 


1  196541.09 
3  196592.88 
5  196712.17 
3  197859.28 
3  202169.9 
3  202714.94 
5  202765.86 
7  202862.06 
1  203164.7 
3  203188.8 
5  203259.7 
3  203532.8 
5  205350.7 
3  205982.1 
5  206038.1 
7  206108.7 
1  206327.5 
5  209675.3 
7  209739.5 
9  209825.3 
5  209926.92 
3  210239.8 
5  210266.3 
7  210301.9 
5  210705.4 
3  210751.5 

I  210777 

7  211030.9 
5  211033.71 
7  211057.07 
9  211061.03 
7  211104.8 
7  211288.02 
9  211295.65 

II  211390.77 
3  211335.5 

9  211402.89 
7  211411.25 
5  211416.2 
3  211487.28 
5  211491.16 
1  211750.2 


113 


3  211780.6 
5  211828.8 

I  214212.4 
3  214258.2 
5  214385.3 
3  214828 
15  220717 
12  221070.2 
9  221074.3 
7  221137.6 
7  221227.7 
9  221232.7 

II  221302.2 
9  221312.1 

1  224027.1 
3  224042.9 
5  224072.3 
7  224115.4 
9  224169.3 
3  225987.1 
5  226011.2 
7  226055.2 
5  230223 
-1  -1 


114 


B.12  N_2.in 


382625.5 
2  0.000 
4  174.5 
2  57192.1 
4  57252 
6  57333.2 
6  101023.8 
4  101031.5 
2  131003.5 
2  145876.1 
4  145986.5 
4  186802.3 
6  203072.3 
4  203088.9 
2  221302.4 
2  230404.5 
4  230408.6 
2  245665.7 
4  245-^01.7 
4  267238.5 
6  267244.4 
2  287535.6 
4  287598.1 
6  287713.9 
2  297150.2 
4  297263.1 
2  301088.2 
2  309132.6 
4  309185.8 
2  309662.8 
4  309698.3 
6  309760.5 
8  309856.7 
2  311691.3 
4  311716.1 
4  314224 
2  317299.9 
4  317343.4 
6  317402.3 
4  317750.8 
6  317781.8 


115 


14  320287.5 
4  320977.4 
6  321065.8 
2  327056.8 
4  330238.4 
6  330273.5 
8  330325.3 
10  330396.7 
2  332796.6 
4  332810 
6  332832 
8  332860.3 
2  333713.1 
4  334542.2 
6  334568.9 
6  336213.4 
4  336268 
2  336303.1 
6  339744.4 
8  339855.7 
4  341946.2 
6  341947.9 
4  342693 
2  342763.7 
14  342752 
18  343116 
10  354517 
14  354955.7 
18  355214 
2  368525.6 
4  368588.3 
6  368704.8 
4  373342 
6  373376 
2  374747.4 
4  374805.3 
2  376756.6 
4  376803.3 
6  376863.8 
8  376953.3 
2  377591 
4  377608 


B.13  N_3.in 


624851 
1  0.000 
1  67136.4 
3  67199.6 
5  67343.8 
3  130695 
1  175463.5 
3  175536.7 
5  175661.5 
5  188885 
1  235370 
3  377206 
1  388858 
3  404521 
1  405893.2 
3  405909 
5  405944 
3  419967.8 
5  419971.3 
7  419979.4 
5  429158 
1  465223 
3  465300.6 
5  465463.4 
3  473032 
3  480880 
5  484394 
7  484525 

3  487542 

4  494240 

5  494338 
5  498315 
5  499708 
21  499851 
9  503625 
3  505487 
5  505518 
7  505561 
7  506292 
3  507022 
15  511384 


118 


5  511440 

4  511493 

5  514638 
5  516631 
7  516639 
9  516650 
3  519414 

7  521868 
3  550218 
15  552731 
27  554419 
15  574940 
5  591043 

8  593665 
7  593704 
-1  -1 


B.14  O_0.in 


109836.7 
5  0. 

3  158.5 

1  226.5 

5  15867.7 

1  33792.4 

5  73767.81 

3  76794.69 

3  86625.35 

5  86627.37 

7  86631.07 

5  88630.84 

3  88630.30 

1  88631.00 

5  95476.43 

3  96225.5 

9  97420.24 

7  97420.37  multiple? 
5  97420.50  mult 
7  97488.14  mult 
3  99092.64 
5  99093.31 
7  99094.52 
5  99680.4 
7  101135.04 
5  101147.21 
3  101155.10 
5  102116.21 
3  102411.65 
5  102661.63 
9  102865.09  multiple? 
7  102908.14  multiple? 
5  103869.4  mult 
5  105019.0 
3  105164.90 
9  105385.3  multiple? 

7  105408.58  mult 
5  105911.3 
5  106545.1 
3  106627.9 
9  106751.2 


120 


7  106765.8 
5  107445.4 
3  107497.1 
9  107573.1 
7  107582.7 
5  108021.4 
3  108057.6 
9  108105.7 
7  108116.6 
5  108412.0 
3  108436.1 
9  108470.2 
7  108477.8 
5  108688.4 
3  108707.3 
9  108731.5 
7  108734.4 
-1  -1 


121 


B.15  0-1. in 

283550.9 
4  0. 

6  26808.4 
4  26829.4 
4  40466.9 
2  40468.4 
6  119837.7 
4  120001.1 
2  120083.5 
6  165987.7 
4  165996.0 
2  185235.36 
4  185340.68 
6  185499.20 
2  188888.38 
4  189068.37 
2  195710.4 
2  203942.21 
2  206730.80 
4  206786.34 
6  206877.90 
8  207002.52 
6  206971.3 
4  206972.3 
2  208346.17 
4  208392.27 
6  208484.24 
4  211521.98 
6  211712.66 
4  212161.94 
4  212593.2 
2  212762.4 
2  214169.74 
4  214229.48 
2  226851. 

6  228723.3 
8  228746.9 
6  229946.6 
4  229968.2 
4  231296.05 
6  231350.08 


122 


8  231427.99 
10  231530.26 
6  232462.83 
4  232536.06 
2  232602.57 
2  232480.1 
4  232526.7 
2  232711.70 
4  232745.98 
6  232747.51 
8  232753.86 
6  232796.27 
8  232959.26 
4  233430.10 
2  233544.09 
4  234402.48 
6  234454.45 
2  238626.32 
4  238731.54 
6  238892.96 
2  240328.75 
4  240516.28 
6  245395.5 

2  245767.80 
4  245816.29 
6  245902.85 

3  246028.95 

4  248009.1 
6  248185.3 
2  248425.35 
4  248514.23 
6  250251. 

8  251220.9 
6  251224.1 
10  252607.7 
8  252608.9 
4  253046.23 
6  253048.35 
2  253789.51 
4  253791.87 
8  254481.5  multiple? 
10  254590.7 


123 


10  254895.2  multiple? 
4  254982.2 
6  255104.6 
4  255140.9 
2  255162.6 
4  255172.5 
2  255281.4 
6  255301.3 
8  255465.2 
2  255622.4 
6  255689.6 
4  255812.2 
8  255691.4 
6  255813.1 
4  255913. 

2  255912.0 
6  255755.8 
8  255759.4 
10  255827.6 
12  255977.5 
8  255829.4 
10  255983.6 
4  255843.1 
6  255897.2 
4  256083.5 
6  256087.6 
8  256123.1 
10  256136.2 
6  256125.8 
8  256143.3 
2  257693.7 
4  257797.9 
6  257963.8 
2  258408.6 
4  258601.7 
6  259286.2 
4  259287.0 
4  260959. 

6  261042. 

8  261180. 

4  261261.7 
6  261354.3 


124 


4  261697.5 
6  261869.4 
10  265220.3 
6  265431.5 
6  265468.2 
8  265578. 

8  265639. 

6  265705. 

4  265762. 

2  265859. 

6  265665. 

8  265691. 

10  265761. 
12  265925. 

8  265763.0 
10  265930.2 
6  265856. 

4  265928. 

6  265961. 

8  265985. 

10  265999. 

6  265988. 

8  265999. 

4  267763.39 
6  267770.85 
8  267783.40 
C  274739.2 
8  274782.4 
10  274920. 

6  275611. 

18  275841.3 
14  275879.6 
10  275951. 

2  275997. 

10  276066.3 
22  276109.1 
6  276263.9 
10  278140 
-1  -1 


125 


B.16  0_2.in 


443193.5 

1  0. 

3  113.4 
5  306.8 
5  20271.0 
1  43183.5 
5  60312.1 
7  120025.4 
5  120052.6 
3  120058.5 
5  142381.7 
3  142382.8 
1  142396.9 
5  187049.4 
3  197086.7 
3  210458.5 
1  267257.29 
3  267375.65 
5  267632.59 
3  273080.07 
5  283758.9 
3  283976.6 
1  284073.3 
3  290956.62 
3  293865.26 
5  294001.60 
7  294221.65 
3  297557.50 
5  298289.4 
1  300228.21 
3  300310.31 
5  300440.85 
5  306584.8 
1  313801.07 
5  324462.46 
7  324658.25 
9  324836.41 
5  324734.22 
3  327227.94 
5  327277.18 
7  327350.90 


126 


5  329467.98 
3  329581.98 
1  329643.43 
7  331820.2 
3  332777.1 
3  338565.87 
5  338690.34 
7  338851.50 
1  343302.6 
1  350026.1 
3  350122.9 
5  350302.3 
1  356732. 

3  356838. 

5  357111. 

3  358667.4 
3  363266.8 
1  365515.76 
3  365550.60 
5  365619.12 
7  365719.16 
9  365846.46 
3  365723.9 
3  366486.91 
5  366594.01 
7  366801.04 
3  367952.20 
3  368526.37 
5  368583.63 
7  368684.75 
1  370326.7 
3  370415.7 
5  370524.2 
5  370900.6 
1  373046.2 
3  374575. 

5  374662.5 
7  374798.6 
5  376067.66 
5  377375. 

5  377687. 

5  378408.5 


3  378420.9 

I  378438.1 
3  379232. 

5  379293. 

7  379356. 

5  380706. 

7  380782. 

3  381086. 

5  392221. 

3  392778. 

3  394090. 

5  394126. 

7  394195. 

3  394516.45 
5  394555.15 
7  394612.70 
9  394688.44 

II  394780.47 
1  398135.0 

3  398131.4 
5  398127.3 
7  398137.4 
9  398218.8 
7  398474.3 
5  398544.3 
3  398582.8 
5  400354.8 
3  400464.7 
1  400518.4 
5  401379. 

7  401475.4 
9  401609.1 
5  401530. 

5  401787. 

7  402530. 

7  403374. 

3  403526. 

3  405805.1 
5  405834.1 
7  405833.0 
5  414675. 

7  415181. 


7  422977. 

7  424998. 

5  426338. 

3  428487. 

5  428606. 

7  428769. 

3  430025. 

3  437015.0 
3  438241.0 
5  438303.2 
7  438395.2 
9  438517.5 
3  439278.1 
5  439329.5 
7  439427.6 
7  442710. 
-1  -1 


129 


B.17  0_3.in 

624396.5 
2  0.0 
4  386.5 
2  71177.0 
4  71308.4 
6  71492.9 
6  126936.3 
4  126950.3 
2  164366.9 
2  180481.3 
4  180724.6 
4  231275.1 
6  255156.7 
4  255186.0 
2  289016.1 
4  289024.0 
2  357614.8 
2  390161.1 
4  390248.2 
4  419533.5 
6  419550.2 
2  438588.5 
4  438723.6 
6  438970.5 
2  452808.0 
4  453073.0 
2  467231.1 
4  467346.5 
2  468075.4 
4  468154.2 
6  468289.7 
8  468499.4 
4  474217.8 
2  478587.7 
4  478682.2 
6  478811.3 
4  482667.5 
6  482923.1 
2  485823.1 
2  492880. 

4  494907.5 


130 


6  494986.3 
8  495098.7 
10  495252.8 
2  499506.4 
4  499535.3 
6  499582.0 
8  499646.6 
4  501511.3 
6  501566.4 
6  503834.5 
4  503947.9 
2  504021.7 
4  510560. 

6  510567. 

6  510746.1 
8  510978.5 
4  514217. 

2  514368. 

2  518684. 

4  518690. 

2  539368. 

4  547311. 

6  547336. 

2  549792. 

4  549855. 

6  552034. 

14  552490 
2  554461 
2  568638. 

4  568773. 

6  569020. 

8  570791. 

2  573696. 

4  573907. 

6  574373. 

2  575204. 

4  575373. 

4  575819. 

6  575853. 

2  576591. 

4  576735. 

6  576947. 


131 


2  581721. 
4  581743. 
4  584552. 
6  584768. 
14  587850 
2  580071. 
8  591767. 
6  592999. 
4  593627. 
6  593708. 
6  594007. 
8  594080. 
4  594337. 
6  594542. 
6  596299.. 
8  596477. 
2  597254,. 
14  597352 
4  597726. 
2  597863. 
4  600092. 
6  600106. 
8  602977. 
6  602434. 
6  615431. 
4  615460. 
4  616588. 
-1  -1 


132 


B.18  0-4. in 

918702. 

1  0. 

1  82121.2 
3  82257.9 
5  82564.1 
3  158798. 

1  213641.7 
3  213797.4 
5  214066.2 
5  231722. 

1  287909. 

3  547150.0 
1  561278. 

3  580826. 

1  582983.6 
3  583019.9 
5  583097.2 
3  600925.5 
5  600936.6 
7  600956.1 
5  612617. 

1  653099.7 
3  653262.2 
5  653605.0 
3  664486. 

3  672695. 

3  677333. 

5  677532. 

7  677847. 

3  684124. 

1  689585.6 
3  689699.6 
5  689890.3 
5  694646. 

5  697170. 

3  704360. 

5  704424. 

7  704527. 

1  707630. 

5  708154. 

3  708296. 


133 


1  708379. 
7  712967 
3  709277. 
3  722666. 
1  731667. 
3  736108. 
5  736126. 
3  737883. 
3  742401. 
5  742407. 
7  742421. 
5  746280. 
7  749857. 
3  796263. 
3  802452. 
7  806625. 
5  808351. 
3  824280. 
3  829588. 
3  831047. 
5  831213. 
7  831504. 
3  832251. 
3  835151. 
5  835321. 
5  837834. 
5  837864. 
3  839616. 
7  840832. 
7  841220. 
3  841280. 
5  841374. 
7  841497. 
5  842105. 
5  843290. 
3  843397. 
1  843449. 
7  847129. 
3  847465. 
3  860874. 
7  861975. 
5  862419. 


134 


3  874447. 
7  875365. 
3  898580. 
7  899671. 
5  901344. 
5  902442. 
5  902592. 
7  904497. 
7  906404. 
-1  -1 


